{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Softmax regression (or multinomial logistic regression) is a generalization of logistic regression to the case where we want to handle multiple classes. \n",
    "In logistic regression we assumed that the labels were binary: $y^{(i)} \\in \\{0,1\\}$. We used such a classifier to distinguish between two kinds of hand-written digits. \n",
    "Softmax regression allows us to handle $y^{(i)} \\in \\{1,\\ldots,K\\}$ where K is the number of classes.Recall that in logistic regression, we had a training\n",
    "set $\\{ (x^{(1)}, y^{(1)}), \\ldots, (x^{(m)}, y^{(m)}) \\}$ of m labeled examples, where the input features are $x^{(i)} \\in \\Re^{n}$. With logistic regression, we were in the binary classification setting, \n",
    "so the labels were $y^{(i)} \\in \\{0,1\\}$. Our hypothesis took the form:\n",
    "\n",
    "\\begin{align}\n",
    "h_\\theta(x) = \\frac{1}{1+\\exp(-\\theta^\\top x)}\n",
    "\\end{align}\n",
    "\n",
    "and the model parameters θ were trained to minimize the cost function\n",
    "\\begin{align}\n",
    "J(\\theta) = -\\left[ \\sum_{i=1}^m y^{(i)} \\log h_\\theta(x^{(i)}) + (1-y^{(i)}) \\log (1-h_\\theta(x^{(i)})) \\right]\n",
    "\\end{align}\n",
    "In the softmax regression setting, we are interested in multi-class classification (as opposed to only binary classification), and so the label y can take on K different values, rather than only two. Thus, in our training set $\\{ (x^{(1)}, y^{(1)}), \\ldots, (x^{(m)}, y^{(m)}) \\}$, we now have that $y^{(i)} \\in \\{1, 2, \\ldots, K\\}$. (Note that our convention will be to index the classes starting from 1, rather than from 0.) For example, in the MNIST digit recognition task, we would have K=10 different classes.\n",
    "\n",
    "Given a test input x, we want our hypothesis to estimate the probability that $P(y=k | x)$ for each value of k=1,…,K. I.e., we want to estimate the probability of the class label taking on each of the K different possible values. Thus, our hypothesis will output a K-dimensional vector (whose elements sum to 1) giving us our K estimated probabilities. Concretely, our hypothesis $h_\\theta(x)$ takes the form:\n",
    "\n",
    "\\begin{align}\n",
    "h_\\theta(x) =\n",
    "\\begin{bmatrix}\n",
    "P(y = 1 | x; \\theta) \\\\\n",
    "P(y = 2 | x; \\theta) \\\\\n",
    "\\vdots \\\\\n",
    "P(y = K | x; \\theta)\n",
    "\\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "=\\frac{1}{ \\sum_{j=1}^{K}{\\exp(\\theta^{(j)\\top} x) }}\n",
    "\\begin{bmatrix}\n",
    "\\exp(\\theta^{(1)\\top} x ) \\\\\n",
    "\\exp(\\theta^{(2)\\top} x ) \\\\\n",
    "\\vdots \\\\\n",
    "\\exp(\\theta^{(K)\\top} x ) \\\\\n",
    "\\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Here $\\theta^{(1)}, \\theta^{(2)}, \\ldots, \\theta^{(K)} \\in \\Re^{n}$ are the parameters of our model. Notice that the term $\\frac{1}{ \\sum_{j=1}^{K}{\\exp(\\theta^{(j)\\top} x) } }$ normalizes the distribution, so that it sums to one.\n",
    "\n",
    "For convenience, we will also write $\\theta$ to denote all the parameters of our model. When you implement softmax regression, it is usually convenient to represent θ as a n-by-K matrix obtained by concatenating $\\theta^{(1)}, \\theta^{(2)}, \\ldots, \\theta^{(K)} $ into columns, so that\n",
    "\n",
    "\\begin{align}\n",
    "\\theta = \\left[\\begin{array}{cccc}| & | & | & | \\\\\n",
    "\\theta^{(1)} & \\theta^{(2)} & \\cdots & \\theta^{(K)} \\\\\n",
    "| & | & | & 1 \\end{array}\\right]\n",
    "\\end{align}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 公式推导\n",
    "\n",
    "Let’s look at one more example of a GLM. Consider a classification problem in which the response variable y can take on any one of k values, so y ∈ {1, 2, . . . , k}. For example, rather than classifying email into the two classes\n",
    "spam or not-spam—which would have been a binary classification problem we might want to classify it into three classes, such as spam, personal mail, and work-related mail. The response variable is still discrete, but can now take on more than two values. We will thus model it as distributed according to a multinomial distribution.\n",
    "\n",
    "Let’s derive a GLM for modelling this type of multinomial data. To do so, we will begin by expressing the multinomial(多项的) as an exponential family distribution.\n",
    "\n",
    "To parameterize a multinomial over k possible outcomes, one could use k parameters $φ_1 , . . . , φ_k$ specifying the probability of each of the outcomes.However, these parameters would be redundant, or more formally, they would not be independent (since knowing any\n",
    " k − 1 of the $φ_i$ ’s uniquely determines the last one, as they must satisfy $\\sum_{i=1}^k\\phi_i=1)$. So, we will instead parameterize the multinomial with only k−1 parameters,$ \\phi_1 , . . . , \\phi_{k−1} $, where $\\phi_i = p(y = i; \\phi)$, and $ p(y = k;\\phi) = 1 −\\sum_{i=1}^{k-1}\\phi_i $. For notational convenience,\n",
    "we will also let $\\phi_k = 1 −\\sum_{i=1}^{k-1}\\phi_i$, but we should keep in mind that this is not a parameter, and that it is fully specified by $\\phi_1 , . . . , \\phi_{k−1}$ .\n",
    "To express the multinomial as an exponential family distribution, we will define $T (y) \\in \\mathbb{R}^{k−1}$ as follows:\n",
    "\n",
    "\\begin{align}\n",
    "T(1)=\\begin{bmatrix}\n",
    "1 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\end{bmatrix},\n",
    "T(2)=\\begin{bmatrix}\n",
    "0 \\\\\n",
    "1 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\end{bmatrix},\n",
    "T(3)=\\begin{bmatrix}\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "1 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\end{bmatrix},\n",
    "\\cdots,\n",
    "T(k-1)=\\begin{bmatrix}\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "1 \\end{bmatrix},\n",
    "T(k)=\\begin{bmatrix}\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Unlike our previous examples, here we do not have T (y) = y; also, T (y) is now a k − 1 dimensional vector, rather than a real number. We will write $(T (y))_i$ to denote the i-th element of the vector T (y).\n",
    "\n",
    "We introduce one more very useful piece of notation. An indicator function 1{·} takes on a value of 1 if its argument is true, and 0 otherwise(1{True} = 1, 1{False} = 0). For example, 1{2 = 3} = 0, and 1{3 =5 − 2} = 1. So, we can also write the relationship between T (y) and y as $(T (y))_i = 1\\begin{Bmatrix}y=i\\end{Bmatrix}$. (Before you continue reading, please make sure you understand why this is true!) Further, we have that $E[(T (y))_i] = P (y = i) = \\phi_i$ .\n",
    "\n",
    "We are now ready to show that the multinomial is a member of the exponential family. We have:\n",
    "\n",
    "\\begin{align}\n",
    "p(y;\\phi)=\\phi_1^{1 \\left\\{ y=1 \\right\\}}\\phi_2^{1 \\left\\{ y=2 \\right\\}}\\cdots\\phi_k^{1 \\left\\{ y=k \\right\\}}\n",
    "=\\phi_1^{1 \\left\\{ y=1 \\right\\}}\\phi_2^{1 \\left\\{ y=2 \\right\\}}\\cdots\\phi_k^{1 -\\sum_{i=1}^{k-1}1\\left\\{y=i\\right\\}}\n",
    "=\\phi_1^{(T(y))_1}\\phi_2^{(T(y))_2}\\cdots\\phi_k^{1 -\\sum_{i=1}^{k-1}(T(y))_i}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "= exp((T(y))_1log(\\phi_1/\\phi_k)+exp((T(y))_1log(\\phi_2/\\phi_k)+\\cdots+exp((T(y))_{k-1}log(\\phi_{k-1}/\\phi_k)+log(\\phi_k))\n",
    "=b(y)exp(\\eta^TT(y)-a(\\eta))\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "\\eta = \\begin{bmatrix}\n",
    "log(\\phi_1/\\phi_k)\\\\\n",
    "log(\\phi_2/\\phi_k)\\\\\n",
    "\\vdots\\\\\n",
    "log(\\phi_{k-1}/\\phi_k) \\end{bmatrix},\n",
    "a(\\eta) = -log(\\phi_k),\n",
    "b(y) = 1\n",
    "\\end{align}\n",
    "\n",
    "This completes our formulation of the multinomial as an exponential family distribution.The link function is given (for i = 1, . . . , k) by\n",
    "\n",
    "$\\eta_i = log\\frac{\\phi_i}{\\phi_k}$\n",
    "For convenience, we have also defined $η_k = log(φ_k/φ_k ) = 0.$ To invert the link function and derive the response function, we therefore have that\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "e^{\\eta_i} = \\frac {\\phi_i} {\\phi_k},\n",
    "\\phi_ke^{\\eta_i} = \\phi_i,\n",
    "\\phi_k\\sum_{i=1}^ke^{\\eta_i}=\\sum_{i=1}^k\\phi_i=1\n",
    "\\end{align}\n",
    "\n",
    "This implies that $\\phi_k=1/\\sum_{i=1}^ke^{\\eta_i}$, which can be substituted back into  to give the response function\n",
    "\n",
    "\\begin{align}\n",
    "\\phi_k=\\frac {e^{\\eta_i}}{\\sum_{i=1}^ke^{\\eta_i}}\n",
    "\\end{align}\n",
    "\n",
    "This function mapping from the $\\eta$’s to the $\\phi$’s is called the softmax function.\n",
    "\n",
    "To complete our model, we use Assumption 3, given earlier, that the $\\eta_i$’s\n",
    "are linearly related to the x’s. So, have $\\eta_i = \\theta_i^T x (for i = 1, . . . , k − 1)$,\n",
    "where $\\theta_1, . . . , \\theta_k−1 \\in \\mathbb{R}^{d+1} $ are the parameters of our model. For notational\n",
    "convenience, we can also define $\\theta_k = 0$, so that $\\eta_k = \\theta_k^T x = 0$, as given\n",
    "previously. Hence, our model assumes that the conditional distribution of y\n",
    "given x is given by\n",
    "\n",
    "\\begin{align}\n",
    "p(y=i|x;\\theta)=\\phi_i = \\frac {e^{\\eta_i}}{\\sum_{j=1}^k e^{\\eta_j}}=\\frac {e^{\\theta_i^Tx}}{\\sum_{j=1}^k e^{\\theta_j^Tx}}\n",
    "\\end{align}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题求解\n",
    "We now describe the cost function that we’ll use for softmax regression. In the equation below, 1{⋅} is the ”‘indicator function,”’ so that 1{a true statement}=1, and 1{a false statement}=0. For example, 1{2+2=4} evaluates to 1; whereas 1{1+1=5} evaluates to 0. Our cost function will be:\n",
    "\n",
    "\\begin{align}\n",
    "J(\\theta) = - \\left[ \\sum_{i=1}^{m} \\sum_{k=1}^{K}  1\\left\\{y^{(i)} = k\\right\\} \\log \\frac{\\exp(\\theta^{(k)\\top} x^{(i)})}{\\sum_{j=1}^K \\exp(\\theta^{(j)\\top} x^{(i)})}\\right]\n",
    "\\end{align}\n",
    "\n",
    "The softmax cost function is similar, except that we now sum over the K different possible values of the class label. Note also that in softmax regression, we have that\n",
    "\n",
    "\\begin{align}\n",
    "P(y^{(i)} = k | x^{(i)} ; \\theta) = \\frac{\\exp(\\theta^{(k)\\top} x^{(i)})}{\\sum_{j=1}^K \\exp(\\theta^{(j)\\top} x^{(i)}) }\n",
    "\\end{align}\n",
    "\n",
    "We cannot solve for the minimum of J(θ) analytically, and thus as usual we’ll resort to an iterative optimization algorithm. Taking derivatives, one can show that the gradient is:\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla_{\\theta^{(k)}} J(\\theta) = - \\sum_{i=1}^{m}{ \\left[ x^{(i)} \\left( 1\\{ y^{(i)} = k\\}  - P(y^{(i)} = k | x^{(i)}; \\theta) \\right) \\right]  }\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后一个公式的推导思路如下：\n",
    "\n",
    "\\begin{align}\n",
    "J(\\theta) = - \\left[ \\sum_{i=1}^{m} \\sum_{k=1}^{K}  1\\left\\{y^{(i)} = k\\right\\} \\log \\frac{\\exp(\\theta^{(k)\\top} x^{(i)})}{\\sum_{j=1}^K \\exp(\\theta^{(j)\\top} x^{(i)})}\\right]\n",
    "\\end{align}\n",
    "\n",
    "以上损失函数对参数向量k求导时可以写成如下的组合：\n",
    "\n",
    "\\begin{aligned}\n",
    "& \\nabla_{\\theta^{(k)}} J(\\theta)=\\nabla_{\\theta^{(k)}}\\left(-\\sum_{i=1}^m 1\\left\\{y^{(i)}=k\\right\\} \\log \\frac{\\exp \\left(\\theta^{(k) T} T^{(i)}\\right)}{\\sum_{j=1}^k \\exp \\left(\\theta^{(j)} T^{(i)} x^{(i)}\\right)}  +\\left(-\\sum_{i=1}^m 1\\left(y^{(i)} \\neq k\\right) \\log \\frac{\\exp \\left(\\theta^{\\left(y^{(i)}\\right) \\top} \\chi^{(i)}\\right.}{\\sum_{j=1}^k \\exp \\left(\\theta^{(j)^{\\top}} \\chi^{(i)}\\right)}\\right)\\right)\n",
    "\\end{aligned}\n",
    "\n",
    "当样本的真值为k时的对参数k向量梯度求导：\n",
    "\n",
    "\\begin{aligned}\n",
    "& \\nabla_{\\theta^{(k)}} \\log \\frac{\\exp \\left(\\theta^{(k) T} x^{(i)}\\right)}{\\sum_{j=1}^k \\exp \\left(\\theta^{(j) T} x^{(i)}\\right)} \\\\\n",
    "= & \\nabla_{\\theta^{(k)}}\\left(\\log \\exp \\left(\\theta^{(k) T} x^{(i)}\\right) -\\log\\left(\\sum_{j=1}^k \\exp \\left(\\theta^{(j) T} x^{(i)}\\right)\\right)\\right. \\\\\n",
    "= & \\nabla_{\\theta^{(k)}}\\left(\\theta^{(k) T} x^{(i)}\\right)-\\nabla_{\\theta^{(k)}} \\log \\left(\\sum_{j=1}^k \\exp \\left(\\theta^{(j) T} x^{(i)}\\right)\\right) \\\\\n",
    "= & x^i-\\frac{\\exp \\left(\\theta^{(k) T} x^{(i)})\\right.}{\\sum_{j=1}^k \\exp \\left(\\theta^{(i) T} x^{(i)}\\right)} x^{(i)}\\\\\n",
    "= & x^i\\left(1-P\\left(y^{(i)}=k \\mid x^{(i)} ; \\theta\\right)\\right)\n",
    "\\end{aligned}\n",
    "\n",
    "当样本的真值为非k时对参数k向量梯度求导：\n",
    "\n",
    "\\begin{aligned}\n",
    "& \\nabla_{\\theta(k)} \\log \\frac{\\exp \\left(\\theta^{\\left(y^{(i)}\\right)^{\\top}} x^{(i)}\\right)}{\\sum_{j=1}^k \\exp \\left(\\theta^{(j 3 T} x^{(i)}\\right)}\\left(y^{(i)} \\neq k\\right) \\\\\n",
    "& =-\\frac{\\theta^{(k)^{\\top} x(i)}}{\\sum_{j=1}^k \\exp \\left(\\theta^{(j)^{\\top} x(i)}\\right)} \\cdot x_i \\\\\n",
    "& =-x_i \\cdot P\\left(y^{(i)}=k \\mid x^{(i)} ; \\theta\\right)\n",
    "\\end{aligned}\n",
    "\n",
    "可以直观的理解为当样本真值分类为k时增大为k类的概率的参数更新，不为k时减小为k类概率的参数更新\n",
    "\n",
    "综合起来的写法：\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla_{\\theta^{(k)}} J(\\theta) = - \\sum_{i=1}^{m}{ \\left[ x^{(i)} \\left( 1\\{ y^{(i)} = k\\}  - P(y^{(i)} = k | x^{(i)}; \\theta) \\right) \\right]  }\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(600, 2) (600,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAASk9JREFUeJzt3Xt4FPW9P/B3NpAlhCQQCLmQEAIpoih4w1QsCkpLqeWRc7F69FioHFBEK4a2GrVQPEgUAWk5HpAWgR6lF3+nxlJbfTheAIsXBK2CggYBIZAAQrIhQALZ+f2xmcns7MzOzO7cdvf9eh4fyGZ2Z7KJzDvf7+f7+aYJgiCAiIiIyAU+ty+AiIiIUheDCBEREbmGQYSIiIhcwyBCRERErmEQISIiItcwiBAREZFrGESIiIjINQwiRERE5Jpubl9ANMFgEIcPH0Z2djbS0tLcvhwiIiIyQBAEtLS0oLi4GD5f9DEPTweRw4cPo7S01O3LICIiohgcPHgQJSUlUY/xdBDJzs4GEPpCcnJyXL4aIiIiMiIQCKC0tFS6j0fj6SAiTsfk5OQwiBARESUYI2UVLFYlIiIi1zCIEBERkWsYRIiIiMg1DCJERETkGgYRIiIicg2DCBEREbmGQYSIiIhcwyBCRERErmEQISIiItcwiBAREZFrPN3inYiIyKuCTXOA83sjP9FtCHy9lzh/QQmKQYSIiCgW5/cC5z91+yoSHqdmiIiIyDUMIkREROQaBhEiIiJyDYMIERERuYZBhIiIiFzDVTNERORJnl8e222IucdJFYMIERF5k8eXx3oiDCUBTs0QERGRaxhEiIiIyDUMIkREROQaBhEiIiJyDYtViYjIm1xelaK5aqfzGlisag0GESIi8iTXb/QeX7WTLDg1Q0RERK7hiAgREZEHpcrUUMwjIps3b8akSZNQXFyMtLQ01NbWhn1eEATMnTsXRUVFyMzMxPjx4/HFF1/Ee71ERESpQZwaUv1PI6AkoJiDSGtrK0aOHIlnnnlG9fOLFi3Cr371K6xcuRLvvfcesrKyMGHCBJw9ezbmiyUiIqLkEvPUzMSJEzFx4kTVzwmCgGXLluHRRx/FTTfdBAD47W9/i4KCAtTW1uLWW2+N9bRERETOiLY6R+Vznt8bx6NsqRHZt28fGhoaMH78eOmx3NxcVFZW4p133mEQISIizzMdHrjKJia2BJGGhgYAQEFBQdjjBQUF0ufUtLW1oa2tTfo4EAjYcXlERETkEZ5aNVNTU4P58+e7fRlERKSD0xBkFVuCSGFhIQCgsbERRUVF0uONjY249NJLNZ9XXV2Nqqoq6eNAIIDS0lI7LpGIiOJhYBqCYSVOJmtUEpUtQaS8vByFhYV4/fXXpeARCATw3nvvYebMmZrP8/v98Pv9dlwSERE5jTUTcUmVsBZzEDl16hTq6uqkj/ft24ePPvoIeXl5GDhwIGbPno0FCxbgG9/4BsrLy/Hzn/8cxcXFmDx5shXXTURECS54sgro+DL8QdloScKNqLi8N06iijmIfPDBBxg3bpz0sTilMmXKFKxduxY/+9nP0NraihkzZqCpqQnf+ta38Oqrr6JHjx7xXzURESW+ji+jj5gk2IiKJ8NRAog5iIwdOxaCIGh+Pi0tDY899hgee+yxWE9BRERESc5Tq2aIiChBpMg0RKrs9+ImBhEiIjLN0A04GcJKgk0PJSIGESIisoVuWPH1duQ6vCLhim8dwiBCRES2CLYsBc4fivxEtxL4squAbpcAaXldjwtnAF9P2XHRR1SkVTeJciPn6IoqBhEiIrJH22b1G2/HRUB2FXw5c6I+XTdcnP8E6DgQxwWSFzCIEBGRK8L6iBgY1Qg7XmiNKYS4Oj2SPjjyMaE15pdLlqkeBhEiInKHXh+ReI9XY3Z6xMI2674+S1UfD57Q7jgeVZJM9TCIEBGlIDt+m5ZeM32w5k03LumDYx4JiZUVIwu673W3irjPkcgYRIiIUpEdv03H85q+Iv1DOsNN8ORMoONI6ME4pjYco/O+6NXKJDsGESIicl9aJgDZ6IGvCEgvhS/34YhDfX1WhH0c89SGRwRbloZWEaUoBhEiIrKGWIwp/um/tuvvwhkg2DmKEa224vxeQGiFL2+F9jGdgoEa+HKqge7fiOOiPUBtiXMKYRAhIiJLKOtCdH/L15pWScsCYKC2ouOYsfMonmvqcbIdgwgREVlCLzgEz/wVOPt/oceEM6HA0e2izmNKIp8XR81JxLX4iuDLW2H7staw89pVtCtKklDFIEJElIpM3MQMr7DRK8rM/B6Q+T2zVxob8VrSy0KBJ3gkVEuSlgmkD4AvZw6CgSVARz2AM6Hi186wYsl55WwKDInUKyQaBhEiohRk6iZmU7+KsAZldowepJfBl79R89Py1SrBY9/uqmGxSufUk+57LZyx9rwJhkGEiIjcodagrNsQQBCseX29WpPO8/l6L5GOtVTHAQRPzIQvb4X2vjvK/XVSEIMIERGZI053qLUsjyIYWAKkpakXlxodPZBTqysR+Yq6RjiCp7WPE0NP+mAgLc34uY3qqAtdTgovz9XDIEJERMbpTHdEE9G4q/tVYWEm2LwQSPMD6UPgy5qsW1sR7ebuy1sRer3Ov+tem11FpR0HQtM+3S6xt3A1gTGIEBGRcWrTHQbrO4Inq0IjIp2jHmrNyuR0N8E79RzQvi2ytqNzukX++o5tEKcVnuwYbUkSDCJERGReLKtDxMJUhKZpfDlztAOCrwjo/o3oUxpn/xy1iDZ4MvRcX5+l1i4FFqmEmGRZyeIkBhEiIopOHjDU6kJ8vTv/0Fsd0tpVFNpRH/pTMyB8GhrpiKe2QhZ89ERts54ku9x6FYMIERFFFTVgpJfBl7fWsWsxxUwxbYq3WXcTgwgREcXOgSWyYf1GAMB/raFVKCwOTQwMIkREFD87py+U/UYMjnSE1YgYJIUeef2H2vmE1lAdiw2CgRppH50w3UqSchkwgwgRESUmee2KryjUvh1Q7y9its26rPeI0yMrvpxqR8/nNgYRIqIkZGalh6FpFZXHAZirw5A3DUsfEP46GudW1dkSXa84Vh4gjLZZV4aOiGmhzusSX8+S9znOYxMdgwgRUTIyM1USw7RKLDfDsGDQ2dzM0Ot0dl2VaOwJo3nzzqiURhlUW63L2qwre52otqGXT8motaMXWtWvw+bvSaJiECEiIscEmxcCHQdDYUI+nSLuiHv2bfh6fKtrJEJoBToOaL9ey1JAEELBxsDNW7fG4vwn0Ytq08vCOrVqTdsET8yMfh6SMIgQEVH8jE6xCF/L9oCRjWx0Ttn4enwr7GOkZQHdLtJ+vbbNnX9RtI9XEWxeCF/uwwie+h/4et0ROe0ihh7l+eT0VglJK4Qyda+HQhhEiIgobkanaqw+Lia+vPhfI4WmTuzGIEJElCrSy0zvmOsFqgWjviJDm9mFPaVz7xlfzxtDf6pMqwSPfTv8AWV9ClmOQYSIKNmllwHdKuDro3HjjmXlipPUNoyLcepDf0qlqz4keOq5yKkho0HOiiXEKYJBhIgoGYk3NwMjB15fDmrp9elNqciWGPt63Qn0ulP/NdPLugJMZ1BRK4o19XWkUGhhECEiSkKG+1IY2eXWA8KmZ9IHd02rdBvStYTWgpt3tOZlqp1a08vgy99o6LWDzQuBc+9HfkL+9YjX4fFwaCUGESKiVGDnLrdOUOvngfAbtqGbt1a79k5R28Kr7eZrZq+d4HEWuKpgECEiooQUPDGzawmwvCdJ98tCy3NbVsOXPS3sObp9P9TChkgeYpSdZbmKJmYMIkRE5H3y6RXx5h88Irv5K0PAHcD5XaHDDO0MrF/8aqplPBnGIEJERJ4X041fXLliZLRCbZWLghRoVDrCGr6eDpVmaUlYgGoGgwgRESUmrRt4ej4AA+3cZaRj5SFDSQo0slDT7SIY6eoqncPrtTguYBAhIiJPMLoLcDBQA19Ote4oiVqNiB5DTdJUlutS7BhEiIhSgdeblgFA8LT25+S73La/F+qAmpYFdB8BX+5jqk8xG0KCLUtle9fI+K/tGjHRWq7r4Pur25QtwTCIEBGlgES4QZlq2d5xIBQKNEJITM4fUq8lkY96KJfrdjaMc/T9TbIVOrYFkY6ODvziF7/A888/j4aGBhQXF2Pq1Kl49NFHkabWrpeIiJJKLL+5G3qOGAw6/9R8jjiSYcdohaxeRBqd6bymaE3RKJJtQeTJJ5/EihUrsG7dOgwfPhwffPABfvSjHyE3Nxc//vGP7TotERHBI8P3sfzmbuA5ETd6ref4ikJ/GP16hTNdfw+rAxkQ+lPexVWu40DEQ2GdYBN0ysQptgWRrVu34qabbsKNN4Z2ORw0aBB+97vf4f33VdrbEhGRtZJs+F4UPPU/QPvW0MoWvaWzYrMzAMHTtaH6j/R+gK801PDszF8BoR1Iy4Uvc1zoNbtdBKQXqW4QGBYm5NM1QmtkGNHoBEuRbAsio0ePxqpVq/D5559j6NCh+Mc//oG3334bS5dqD1m1tbWhra1N+jgQCNh1eUREqUX8Dd+DqzyCrc/Dl/XvBo9uN1dLIjq/F2j7S+jv3S4Cet0BnNkQ9lpmplSUxwaPfVt1ZIT02RZEHnroIQQCAQwbNgzp6eno6OjA448/jttvv13zOTU1NZg/f75dl0RElJpMbMwWq7CpILN1Eue/Mn6sIDufIMRXj9HZL0SvxsRYZ9as2K/DrCTbmde2IPLHP/4RL7zwAtavX4/hw4fjo48+wuzZs1FcXIwpU6aoPqe6uhpVVV3NXgKBAEpLS+26RCKi1KC3MZsVNQzxTAUFjxs/VmjqOp8Z3So6m48hclRI69rN7COTPljayTjiHBaPQiVbvYltQeSnP/0pHnroIdx6660AgEsuuQQHDhxATU2NZhDx+/3w+/12XRIRUWpzum7Eqd/cDZzH1/MmoOdN1p5XRjkyw5UzxtkWRE6fPg2fzxf2WHp6OoLBoF2nJCIihwVblka2UhdaAcS6P4xOqPD1D/2ZPhgQThg6j9iJNRhYYnxfmBh5YrVSgrEtiEyaNAmPP/44Bg4ciOHDh+PDDz/E0qVLceedd9p1SiIiEvmvlfXb6Fx+6itC5C61cTp/KPKxjgPmemvIls3q3ax9ve4I/dn5esETU4FgU9ceMZ0raYKBJUBHPYAzADr3jumoD38xAxvdmZakq5XsZFsQWb58OX7+85/jnnvuwdGjR1FcXIy77roLc+fOteuURETUSW3DN1/eCutXd4ghR60OQm1ZqxpfTwSb5sLXW79LavDEzNCyXHGEIdgECK0RK2lURz58/cI/5AZ0nmBbEMnOzsayZcuwbNkyu05BREQKulMDFq/uEG/4WiMewRMzAV/P6K8hGwVRvf7OIlBfdlXXzrhiY7GeUwF0134u0PW1dxuofgHKEKUMUDbuI2N0o79kxr1miIiSid7UgNrIRRw3U72bv+meH8rrTy+L2uvD13Oy7GJOR/3atXqVqIWo4ImZQHp+6PN2hgFO5TCIEBGlEstXc8RxIzXUe0Rv6TEgG+3JjO3cKq8lDz/B5oXAOZWu4NxXxhIMIkREZDvV0GEmxFg9cqDzemKxq6/P0lAISfFRCzsxiBARkf3UNouzYxWPVTrqwzatU9X5eNgSZrO9U6I1O+vctC/ZMYgQEVEYO3phqE1h2LKKx0pG+6Gc+0L6q9n3h1M7DCJERKRk4TRIsGVpaNdbJWklyyXhm/EpRw7MtEfvVgJ0XKR6LsDk5nrpA4COLxE8ORNIrwj1JlGrFRFa41qJZCT0qTaNSyIMIkREycTpDdH0znf+UPSVLMpVMHGsUPFlVwHRbtgmNtcLLUsO70Xiy31YfQSnW2T4MTyqZCT0yUZckhGDCBFREnG674Tu+WRdU9UET1Z11WL0+C58ve4OfzxjjKm27MGToSDi67M0/O8nZoavqtGr+5BfV+fjpvqwqNXEAKElxoDhdvPBk50N3JIYgwgREdlH7yba8WXXiMBZAJ1BRHpca8pGTv45eXgAgLS0ruuQTfPoBqjzn8RVu6Jb+6FsN695XHKHEIBBhIiIvKxzRMXwSE/6YCl8hIUBoVV3dEYUPDFVO4QYaAjHje/MYRAhIqJwTteZRNNRJ/01YrpEzlcEX96KsPAh7UsjtmwXm6OJG+J1KwkrAg09Xhfav0ZZ99EZQKKNdASbF8KX+zC7pZrEIEJERGEs/a091lAj+7x0g5dP40Qc3/XX4MmqqFMrWrUZhmo2Tq0Gzu2KDDEnYqzl8FLocwmDCBER2Ua/mLXV+POM1omkpQHpFaGlwWHHlADoHPlAZPAw1Eb+3C6g7S+hZcLyFTqmQ4jBKSdfke6mgYmOQYSIiBwVbFkaWtYrnAlNl4jTIOKKFfnogmw/F6MjNXrHaY58ODml0nEEwRMzdTcFNL1pYAJiECEiIlPi3bpetzlXR52hFSuh5bln4OujfrM2cp3B5oWArwC+7Gm65wsxVvCqKwVGOoxiECEiInMsGjmQ99KQClHFwlIjOmtGgse+Hd7fw8ymesHjQJrf+EUbWU5roO7DzZGOeIOk1RhEiIg8yms3DMvJe2lEK0TVYqb9uwbL9nrpNkRqVub574vHVvUwiBAReZXHbhheY0WICJ6sAtIHGO/eqtFgzfPhw8MYRIiIyH0ml7FGHS3yX2t8k7iOLyHVfRhYlcPdcq3HIEJERO7oXE4LxDCiEG20yOyUzfk6U9cQtupHvmw3GabLXMAgQkRE5hjt56HD9q3tDe9Pc6Cr4LWz0FXqviodL2tg1raZU2YWYhAhIiJTPP9bv5n9acSGasqVOu1bwsOGsoFZIrMoSFqFQYSIiJKLgS6nem3gk5nXgiSDCBGRV3nsN1dP0aoD0WgZHyEtLbyrqyjV31cXMIgQEXmU135ztYrmihcTxZ5Rd8EN1Og/P8Y9cMh6DCJEROSsOPuj6AUZX051bK8rFah27oGTXhb6RGcRq8R/bdfH8pUzFo2mWBHUEgmDCBERJRadIBMMLAktyTW5tNZoUzPbV/ukWCM7BhGiFFPz77/EV5/VRzw+8MIBqH7+fheuiMhiyhUvURgZfUi1EQqnMYgQpZivPqtH3Yf73L4MIm8wMvqQYiMUTvO5fQFERESUuhhEiIiIyDWcmiEiImeZ3OAu5aTY+8MgQkREMYulkDPuAk8zN+r0ssjltx6XagWwDCJEKWbghQNMPU4UlQuFnLo3al9R6O7mK4Ivb0X0Y42EmhQboXBamiAIgtsXoSUQCCA3NxfNzc3Iyclx+3KIiEgheHyyehDpdhF8/WqtP9+Jmep7yXTumqv6HLVRGyMhhWJm5v7NEREiIkocvp5AUOXxtDTt56iO2nyK4ImpQLdLDDcyI3swiBARUcKwrH4ivQy+vLXWvBbFhUGEiIhST1oWAHZN9QIGESIicpxmAACcDQHsmuo6BhGiJOT2fjJun58cFOuKEgYA6sQgQpSErN5PRgwWpcMG4OEX7tcMGqMmXoo7F9zG/WxSSMwjF+mDAaEV6Dhg7QVRwmEQISJdymChFTRKh7EXSTJwYtpEXGobPPZt+8OIaqOzxGlwluxsDSL19fV48MEH8be//Q2nT59GRUUF1qxZgyuvvNLO0xJRguAUjkdZOG2iWwzaWTRqJxadepttQeTkyZO45pprMG7cOPztb39Dfn4+vvjiC/Tp08euUxJRguEUTgpgLQjpsC2IPPnkkygtLcWaNWukx8rLy+06HREZkF/az+1LCFM6bADOtJxBfV2D25fiaUm3xDS9rGtqJGNM5DRJtxLnrsXi9u1h3yt2bzXEtiDy5z//GRMmTMDNN9+MTZs2YcCAAbjnnnswffp0zee0tbWhra1N+jgQCNh1eURJTWvfmBHXXejwlUT38Auh6ZepQ+9jGIkmmUYV0svgy98ofeh2V1PLg5z8e9V5h/XMUmWPsi2IfPnll1ixYgWqqqrw8MMPY9u2bfjxj3+MjIwMTJkyRfU5NTU1mD9/vl2XRJQy9OorzNZmiMHGaDGq3sZ6C2//JQ7urpfOl5mdaeh1KQnoNRIDku/mnExB0ga2BZFgMIgrr7wSCxcuBABcdtll2LlzJ1auXKkZRKqrq1FVVSV9HAgEUFpaatclEqUss7UZynAyauKlEaHkbGsb8gp7qx6v9Pm2Oo6AeFm0aQmrdpw1cHMOtiwF2jarXoNng4q48y/QNeUkn3rikuUItgWRoqIiXHTRRWGPXXjhhfjf//1fzef4/X74/X67LomILHLngttMHb/6kfVo3H8MANDWehaZ2ZmouKxcqllRhhqtERVyhqU3eWVwMbNs9vyhhBtJUKsJUe4K7MiS5QRiWxC55pprsGfPnrDHPv/8c5SVldl1SqKUozXFAtizBDbW5bbTHo8eXMRaEUo+nh25sJEXliwnEtuCyAMPPIDRo0dj4cKF+MEPfoD3338fq1atwqpVq+w6JVHKcXr5ayznY68QC1i8soNsxpoQU2wLIqNGjcJLL72E6upqPPbYYygvL8eyZctw++2323VKIvIg9gqJXyqOKlDqsLWz6ve//318//vft/MURBQDvVUtK6rW4eNNuzhqQWQHeZ0MR7W41wxRohtQUai6/DXaUlu9cPHeXz7gqhaylxMrczxKWbya6hhEiBxgV1Fpfmk/PFb7s3guDUBXXw8A7HRKjjA03ZSeD3S7KPLxJA8qqYZBhMgBRuoknnt0vellsYMuDrXCjrcg9OBuY3UcelM6crF8PURyvpxqty8hNiwuNoVBhMgjGvYd0z1GDBziSIh4o3eqINTMyI349Rjtsip/nHUplMhYXGwOgwhRAhEDR92H+zB16H0YOqrC8h4c8poToy3d1ZxtDe0bxS6rRBQNgwhRgqqva7BsjxZxdMJozcmKqnU42diEvMLeuHuJ+pYNxw4ex9Sh94VdY+mwAXj4hfulURDWoxARgwgRRYxa/L+nN+BfH5ikefzMperhQ0krZBitSSGi5McgQkQRxBDidAt5KyXytROlEgYRIgdE28RN/Fxheb7u6yhrNsSPzaxmMSORu6Im8rUTpRIGESIHGPnt28hSV63CVKt+uxdHEcRajniohaB4il+JKDkxiBAlCK2phpFjh2sWjJolCKE/xRUv8YgWjoyMEBFRamAQIUoQTkw1WL0UWAvrM4hIxCBClIK8UsipbGYmsqp1PRF5H4MIUQoyOrpi9xSKnct4Of1DlBgYRIhIkxUjI/LRl/zSfuiR5UfBoHxMe/w21eLVMy1nor6GktYIDqd/iBIDgwiRS7RurqMmXoo7F9xm6QqWWMyb/CSOHvw67LFYrkUcfRlQURgx3aL1WivnrFN9DSJKPgwiRC4RV6iIxHbn4ihBPDdfvZBjxNGDX1t68xdbvRsZ3bBqFRAReR+DCJFL1EYDpg69T/N4M03LtEKMF/p4cHSDiOQYRIgcpDUaII4EqG1iJ9ZMWFnzwEJOIvIKBhEiB4nTMWZ2na2vawjbxVas03j5mb9h19bPAYQakB07eFy3hkNsVGYk1KgVjRIRWY1BhMhB8pAwdeh9psKI0q6tn+PN371t6vzHDh4HENm/Qwwwtc/8DZNnTcTKOeuQmZ2JisvKw57v5dES5ddkZcdZIrIPgwiRTWKZhikdNkDa/M6ujezUpKWF/jyy9ygAGLqB63198o+B+OpTjEwlKXuSiOczep1E5A4GESKbxFKUKR8x0btJRtsPRi/ERJu+mTt5EXpk+dHWejZs+a5y2sfo16f8OoyEipef+RtumjVR8zXUaE0lWVUcy0BDZA8GEaIEJU6zyMVb2DpzafhIiJnpI6OMXNvhzpEZPasfWY/G/cfQ1npWmkoyU39jBlf7ENmDQYQoQeWX9lN9/MWlG3Bz1SSsqFqHjzftUn2e2FhM77f8oaMqpCkkK6eE5AFCbdTlZGOTodeZ9rh6T5Roy6CJyFsYRIhcIq+ZiOUmr7cp3MylU6QRjQEVhVKg6JHll44RBKB/aV8ISAsbYRFX99jV0VUeIGIZdYml/oaIvIlBhMglVtzk9W7IY26+Gu1nzkVMuVh5DXrsCA160yTyol8i8jYGEaIEpndDlo886AUCrSW9yueJjxtd1eNGbYU8YDm5+oiIzGMQIbKJF26Azz26Hg37juHhF+7XDQTK5a8i5fOs7vSqNkVl5WiGVdfphe8nUTJiECGySaw3wOceXY87F9wmhQg1Q68cjH99YFLE42ItiHhzl2+sJz4W76oSsdPr0FEVlkztqL2G0Y35nMQlukT2YBAh8pgePXsAMH8zHlBRiLWfLw97TD41o9fVVT4yodfevb6uwXBth1UBKBZGdvolIncxiBB5zG0P/zMA8w20xGBg5OarFiKUIxNzJy8yfe1qYm1rH42X61OIyBwGESKHGf0t3chNdNTES6URh4JB+Yafp0YsVhWvQb7M1wgjX5fYlyS/tB8GXVxi+rXFHihujmRwlIXIWgwiRA6z8rf0WGsp1ApEtYpVrRx9iLWmRBBCUzvv/Hmb6k7ETuIoC5G1GESIEthzj67H/p2HcOzgcVM3ZbXjlHUhfQp6AzBXpClvnCantuGdcpTDyPXa0XKeiNzFIEKUwMQRkalD7wvrAWKEvG+IWiGpsgnavMlPSu3Y1UKPkUAhpxxZiLW2hYgSG4MIURzc3JF14e2/RFpaaMSi7KKSsD1bjNCaipFb/ch6fPDqR4ZWvIg1JUbfE3GURG2pMWDNKhsjO/0SkbsYRIji4Ga9gHwEZH7tg5g3+UkAxm++atMlyumZtLTQn+LOtmqvo2T0PVGOqKhNF8U7FcPCUSLvYxChlODmyIXaOY18zmwnzz5FeQCM33y16klWzlkn/f3OBbcZLog12w01UadiOMpCZC0GEUoJXlrpYDQo6B13puVM2I169orp+OXMVWgNnEHBoPywZmZyK+eswz/e2hXxuHjzv3tJV23IhhWvYdLMCYau1+wKHiPfk/6lfcOOUY7iuHHj5ygLkbUYRIgSgHJDOrF+Qj5dAgD3r5ih+1r/eGuX4VAmhhA7RpTyS/tpfk4MHPNrH8TUofdJjzu9VFckX+Ej769SWJ7vyXb0RImEQYQoARgpLLVz+smOESWjK2wqv38lmo42WXpus8SvX/keVFxWziBCFCfHgsgTTzyB6upq3H///Vi2bJlTpyWylRU7shoJEEbO46XpJ6OM1IkolxETUXJxJIhs27YNzz77LEaMGOHE6YgcY0W9gJEAYcV5Vj+yXrNuxArirsFmwlkihicispbtQeTUqVO4/fbb8etf/xoLFiyw+3REqqwYuXDLe69sR+WNV0TUiYj6l/bF/NoHdV+ncf8xy69NXufRdvocABZzEpE5tgeRWbNm4cYbb8T48eN1g0hbWxva2tqkjwOBgN2XRzF44LVXUHfiRMTjFXl5eHrCjS5ckb5Evjl+tfswKm+8wlCdiBFWhbIBFYVhdR6cQiGiWNgaRH7/+99jx44d2LZtm6Hja2pqMH/+fDsviSxQd+IEdh076vZl2MpLfUcCJ1qkc6tRa0ympq31LAD9UDZ38iIcO3hct0ZFPC93oyWieNgWRA4ePIj7778fGzduRI8ePQw9p7q6GlVVVdLHgUAApaWldl0ikSYv1S60n7FmyuPowa8xd/KiiNUqenvO6J3b7vdKbDOv5GTISeSpPSKvsy2IbN++HUePHsXll18uPdbR0YHNmzfjv/7rv9DW1ob09PSw5/j9fvj9fuVLEaW0k41NACJ7iQDhm88ZuVkeO3hcei0g1JdDb8pnRdU6fLwpsgHa9bePwc1VkzRHZJTt4vWuTetzjfuPuR4KOapDZB/bgsgNN9yATz75JOyxH/3oRxg2bBgefPDBiBBClKri+W07r7C39He9m2V+aT9k9goFfTM79Z5sbIoIAgMqCnFz1SQA0ZuMzZ28KOprG7nBNx9tNnCVRJSobAsi2dnZuPjii8Mey8rKQt++fSMeJ7Kbl+sY9M4tho14u4rKp2TOtJxB2UUlMb+W2Fper5ZG3oU0Fstm/ho7Xv9E/0AiSljsrEqmVeTlmXrcC7xU82GWfO8XNWZD1rKZv0ZBWb6hJb9AqI25spW8OB2j974WDDK3EZ7Snve/iOv5ROR9jgaRt956y8nTkU28ukTXSolUnGg2ZM1eMT3sY72v1cwOvErKBmpqK3KIKLVxRIRIRbLdJOUFp8oQYORrlRe3yv9u5HkHd9dL5zx28Lip0GSkmJWIEhuDCFEKUCs4BbpGKJTkq3GA8OJWM4Wu8TZhS7ZASESRGETINVodWgFvd2m1mpvN08yOUACh6xIEmy6IiFIOgwi5xskOrV4e4k+EQtqBFw7AmVOh7Rfk4SjeWhpxozwiSl0MIpQSEnmIX9nIzI1CT63z6V2HXlOzhn3Wb8RHRImFQYTI4/TqLLwy2qMMTFrt4omI5BhEiBKckdGRPgW9pV4g8nAS69SKvK5FLGw1Upg64rrhAICRY0N/xttnhIgSH4MIUYJT24MGCF/5MnOpelO0WKd45HUt/Uv7AjA2MqO8DmWfESJKPQwi5JponVi93KXVavEWfBoZiVhRtU41jDz36HrVOo3C8nzDRaQC0gAYDzWrH1mPaY/fFjaqkl/aL2JXYCJKDQwi5JpUWZ6rJ96CT70gM2/ykzh68GvVY/TCxv97egP+9YFJUY+R9yHRGp0BIsOG2mohN5cyE5E7GESIPGrh7b/E59vqdAs+9W7QBz49JG1Sp3WettazqnvP6IUQJSOjM9GmYxJhKTMRWYtBhJKK3b9RO/kbe2F5qJBz6KgK1cfVKKdazra2ITM7M+o0j3wEw8u7FBNRcmIQoaRi92/UTv7Grjdtsmzmr3Gm5XRYOIi3OdhXn9XjTMsZ1REUeTdVcdUL0LUTLxFRLBhEiDzIyMjEnve/sPy8RotG716ivgqHiMgsBhEiD3KrVqJocAEA9SCUX9oPPbL8UVfUeKW5GhElDgYRohQn38ROXOKrDEIDKgoNjZQYqSERlxLLg4k4vaO3AsgL7e6JyFoMIkQWMruJm3LkQd6EzClGbuRizUi8xbpzJy/CV58exMylU1SP13sNI6uIiCixMIhQUom3OVi8r292E7dEW64a7/UeO3gc9XUNmDr0voiC2BHXDdfsAGt0KTPAXiREiYZBhJKK3Tcavdc/29oW82sPqCiUpijUVqIkwyZysX5dB3fXG/7aEy3cEaU6BhGiGKh1ENVa9mrEgIpCrP18ufSx1vTM3MmLkNnLDyB0U09Li+l0rtH6uqYOvc/hKyEir2AQIYqBVgdRcYdbs4zWYMgLRs3WkpiZslBOQcXbK0Tv3LEGOCJKfAwiRB5i57SCmde2eopL79x6ozt21/4QkXsYRIhiYNWNUTzere6k8roU0co563CioSnssYJB+Zj2+G22BQK90R0WmRIlLwYRohhYdWM0+jr5pf0sOZ+csi5FFK1rqluBYPUj69G4/xjOtraF7fbrxnJnIrIWgwiRxxltJqYnv7Rf2PSIWl2K2D1VJHZRXXj7L3G2tQ1lw0ui7p5rF/k5pw69L+oKGk7jECUWBhEij9MrZAWM9ch4rPZnqjdxef2Gso6j4rJy3LngNqk499jB47YEEXkjuHgLWzmNQ5RYGESIEoSRYlO1ZcXykDJ0VIV0IzdTl6LXgl2PmUZw7ANClFoYRIg8ZNTESyMCQsGgfMPP11pWLIq1niLeOgw7G8ERUWJjEEkwD7z2CupOnIh4vCIvD09PuNGFKyIrmdmnJtGojdYA8TWCI6LExyCSYOpOnMCuY0fdvgyyiVp9RLKsDIk2WmOmEZx8xIgFqESJj0GEyEPcqI9QrqbxumQIZUTUhUGEKEFE++3f6MjA6kfW44NXP5KeU/38/Rh0cUlYb45YXtcIK66fiJIPgwiRx51pOQPA2LJU8VgtjfuPRYx+3LngNkdqU4wuq2UfEKLUwiBC5HH1dQ2YOvQ+DB1VoTkt8dyj67F/5yFkZmdG1Fsk2g2cfUCIUguDCFECqK9riLqyxOiIBpfJEpHXMIgkmG+WlKg+XpGX5/CVkB3iraMQ92RR6lPQGzOXTtGsBbGTXqdUeVdVIko9DCIJ5pEx49y+BLJRrNMS4s386IFjqr060tLivbLY6a0EkndVJaLUwyDiMVoNywA2LSNt+3ceAqAfZPSKWYmInMYg4jHyhmWDcnsjKyND+pwguHVV5HXHDh7H1KH3RdSRiM3QVj+yHgd2hRezJloRKxElJwYRjxnSJw+t7e0AgDemTHP5aiiRKHfVlbNjx1wiIiswiHjMsu+Gpl7u2lALgHvLEBFRcmMQcYlewMjs3h0A95ahxKfXoKyw3PjuwkSUfGwNIjU1NfjTn/6E3bt3IzMzE6NHj8aTTz6JCy64wM7TJgQGDLKSl7uR6hXQcukuUWqzNYhs2rQJs2bNwqhRo3D+/Hk8/PDD+M53voNPP/0UWVlZdp464ZXk5Lh9CZRA2I2UiBKVrUHk1VdfDft47dq16N+/P7Zv345rr73WzlMnvJ+MHuP2JVCK0mtARkRkJUdrRJqbmwEAeRpdQNva2tDW1tWCOhAIOHJdVhWELtm6BQdVrrk0JwdzTAaL3/5jB3448nJTz3Ebe6AkDrWwIS711WtARkRkJceCSDAYxOzZs3HNNdfg4osvVj2mpqYG8+fPd+qSJFbVa5gNG9HsOHIEPxxp2cuZFks403ofB+X2Zg8UFW6OPDBsEJFXOBZEZs2ahZ07d+Ltt9/WPKa6uhpVVVXSx4FAAKWlpU5cniX0bt6Lt27BoUAAy757o+beMMrHjR5nNavC2aDc3uyHooFhgIjIoSBy77334i9/+Qs2b96MEo1N2wDA7/fD7/c7cUm20Lt5HwoE8HFjqOmU3jTF16dPGzrODDd6koidYTltQ0REamwNIoIg4L777sNLL72Et956C+Xl5XaeLiHsb27C9etWSzfoIX3ypCZmADD71VfwcWNDxHFqx5rl5pJhLle2DotJiSiZ2BpEZs2ahfXr1+Pll19GdnY2GhpCowG5ubnIzMzUeXby2t/cpPm5vSdPSJ+PdlwyKeqV7fYlJBS7p3S83JOEiJKPrUFkxYoVAICxY8eGPb5mzRpMnTrVzlOb4lQdRklODobn9w97bEgfYzUhdlyPFwzK7Y1Vkya7fRmErp15OapCRE6yfWomEThVn/CT0WN0+4N4pVYilnA2tmxQRLDSa8yWyjUkbo48aJ1j5Zx1uHvJFNvPT0Qk4l4zDpr96ivYezL8ZuvVm2ws1xRt+XJRr+yoNSKpWEPi5sgDRz2IyCsYRCykN4qw96S7N1u7pqCMjGZUDhiAI6daIj6vHEEhIqLUwiBiIb1RhNb2doeuRJ3e9c3YUBsWFoyO1ihn4Frb2yMKbaddPgrTLh9l+FqJiCg1MIg4ZMaG2qirYJZs3WKoM6sdvUDkS4Zjobak+Pp1q1Nm1Y/TuKqFiJIJg4hD1KYl5NT2qFGjVUsRzxJY+ZJho/QCkbz/CVmL9R1ElEwYRFwmBgi91SXROLEEVgweYlM1q4tLU23ZMhERhTCIqHCqFbo8QOgt641GbwmsFdctBg8xOGkVmcY6MuPFlUNERGQ/BhEVdiwlVfutXryZ6wWI2a++AkC9FkMu1uvWW1orkgeneFrNExERiRhEYmC0sFQu2m/8Wn3fxMeVvUdipexjIgadC/r2xZFTLSjqlY2e3btLny9VTBeZGXkRQxanVYiIKBoGkRgYLSxVI4aBol7ZhkcXRhYU4h+du/ZGG1nRo9XHZM7oMaaClZGRF46YEBGREQwiDktLC/0pjjw8vuVNPDJmXNTnLLj+27jjTy8CcKaW4tE3NuJUZ8+Tkpwc/GT0GNONx8TAFe+OwURElNwYRBymDBIFWb0AhKZ7xJGW0+fOSct9xemOprazuH7damRlZFh+c//Njm1oOnsWBVlZuGPk5Vhw/bcjjjF7vnj6kpA31fz7L/HVZ/URjw+8cACXFBNRzBhEXPYfnd1GlVMjyoZgQ/rkWVYronUNehZsfhOPXht99EY8LisjA8Pz+7NGJIl89Vk96j7c5/ZlEFGSYRBRYfWeLPLiTnE0w2xDMCMjEnrXPaRPHop6ZSMNwOFTLWHTLUb2izl++rTuNQAwFFaMnpOIiJIbg4gKvRtgflZPU6+nVtxpdqmtfMWL1k1a77qjhRk3dr9NxR13iYgoHINIDPSKS5XkTb5i3W3W7M698poTIFR30rN794jRmFjqTezaxZeIiFIPg4gNlO3Q7Wq/ruwLAnQFi7cO7I8ILsPz+wOIfSTi9LlzANgFlYiIrMMgYpEHXnsFggDNfViU4cQKZkdJ4nXkVAtmbKi1fV8bIiJKHQwiFtEqupR/3kxo8GJn0iF98qQ+KJR6Bl44wNTjRERGMIgoeGUlh3LUpLWzwVg0g3J7SwFGXositm7v17NnxOfM1Kwor0mrZbxRRnfcdWoTQoqOvUKIyA4MIgpWruSIFh70Cj7lG90t2PymanMwZdiQT5lEm/5R+5yRUKDcfC/eqSGjIYKra4iIkheDiI32Nzfhrg21eFYWEMRwoncT/rixQeol0qNbN6nQFOgKBmqBQjl6oNe35JslJXhkzDhDoUB+TURERFZgELGBfHQhU7abLRAKJ2KrdqArKMinOVrb27G/uUkKHz8ZPQY/UdmU7r+3vYvPv/4aQNeeMPKdfFvb26XXNDOqsGDzm1LzstKcHMwZPQYzNtSGXZOXcSqHiChxMIioGJTbWwoKxb2ypTBRkpOj+ZyKvDy0tmsvb9Wa+hCLP2OZ5rhn1DcjHlOOkty1odbQay3eugWbDuyXQpBoeH5/zBk9Bj27d49o2S6fGjJSw+IUTuUQESUOBhEFZa2FUXq/aet9PtYbud5v/8oRGT3iHjEiMWyoTQMpH5thMPQQERGJGEQULujbF4Cx1TNLtm6J2KxOzYwNtdJuunLXlQ3CT0aPweNb3owIAOJ51Lz02S7804XDAVj327/W9I/oNzu2YefR0Hnys3rikTHjpOkk8f0Q3zursZMrEVHyYhBREIOFkRu8vIV6NEdOtai+ltj63WzLeDGEWGnJ1i1468B+aamvSKwRUduh9+PGBgCQ6lKMhLJYsK6DiCh5MYi4qO7E14aOW7D5TTx67biY9oiJVtcidzAQQGt7u+FpqTv+9CIA4I0p0wwdT0REpIZBxAFjywZpNg77723vhhWdKpuEtba3S4WzsUzDiNMtRqY3xPMYmZbq27Mn+nY2SPNKEzj5Oc08TkRE7mEQcYCZKQu11TPxLJn9zY5t+I/LR+mGAXFDO8BY4JGPxnhtlQqncoiIEgeDiAPs7Guh99v/zqNHw/qWyMmnd46capFqVrRGb+xeouu1kRUiIrIfg4gDrBoxkAeE4s7QYGSkQ609vJJ82XK02hNxia7RFUNmeG1khYiI7McgosHI3iulBgtB4yWOVBgpTlVSWzasRm/ZsjgiIa6oMbpiiIiIKJqUCiJmhv6NTAOIIwIzNtSiZ/fuEXu6GF3ZsmTrFhwMBDSPFW/+VvY20fpaOCpBREROSqkgYvVNdvHWLfj866/D6itiOcfBQEBaKVORlxfRy+Pi/qFiVbGFvBqxl4dypEJtZEd8fXFEZ8HmN9GjW7eoDc1EaqNARkaPiIiI1KRUELGakRu3UUZ25TXS46Nf55JakZGRnUevNd5QTRw5ka+yYREpERHFikHEoGjTIt8sKZG6o8oLSisHlBh+fXFX3hEFhVGnePRqOHp00/6W1n62C5Pj7Mo6+9VX8HFjg7QKRz4ttXjrFhzqHJHp17OnqYADcGSFiCgVMYgYJJ9yUe7OK2/RrlbnYbTB1v7mJt3mZfFML8lDiDLQGK1n2XvyBPY3N0m9TeTFsGojRGaWLnNkhYgo9TCImDQot7fptuZmmonFcj1ZGRnSSIxWDxDR//xjB+4YeblmoNELTeLrqx0njoiI+9MALH4lIqLoUiqIWDH0b6YNOqC98658BKKsd640wqAXJOTUQpHehniNra2qjxupUQEiR3yKemWjrHcuAPURkSF98tDa3m6olwkREaWelAoiVg79G/1NX2vnXTm93Xe1OprqhSInupHqFdCKweX6dasZRoiIKILtQeSZZ57BU089hYaGBowcORLLly/HVVddZfdpLSeOVIjLdK0i9v0QN7tTq9XY39yEuzbU4lmNm74V0x/xhhq956m1mAdCxa9paawPISJKVbYGkT/84Q+oqqrCypUrUVlZiWXLlmHChAnYs2cP+vePfSM3N8TS1dQIse+H2mZ3cocNdkiNV6yhJtbnyXcaJiKi1GNrEFm6dCmmT5+OH/3oRwCAlStX4pVXXsFzzz2Hhx56yM5Tx8WNqQ5lsamyZkXtcTP1JEpar09EROQk24JIe3s7tm/fjurqaukxn8+H8ePH45133lF9TltbG9ra2qSPAy7tZ2LlSg+jS3eVIy5agSfeICR2bHV6KkQZdNgXhIiIABuDyPHjx9HR0YGCgoKwxwsKCrB7927V59TU1GD+/Pl2XZIrzCzdFetEinplIw1AZmejMHHJ7eod23Dy7FkAwNnz5w23ZZe7Z9Q3pb//9h87cLRzFU3/rCz8cOTlUmAwu9JFbURFXmSrNbXV2t6uWT9CRETJz1OrZqqrq1FVVSV9HAgEUFpa6uIVaZP/Rq/cG0ZrV967NtSG1Xoob8JinYi8cRoA3DHycgDAtMtH6V5LtMf/e9u7KOqVLS3x/WHn68rJA4OZlS5aQePRNzZKfxeDlkgMO+LSZTm96bFYN/cjIiJvsS2I9OvXD+np6WhsbAx7vLGxEYWFharP8fv98Pv9dl2SpcxObTz6xkZs3Lc34nG1m7BIbPuelZGBygElmi3TjV6LfDREy4wNtcjK6I6nJ9yIEQWFUlDSCjt6gWHB9d+WHktLCz8mKyMDw/P7q7623vSYcnM/IiJKTLYFkYyMDFxxxRV4/fXXMXnyZABAMBjE66+/jnvvvdeu0zpCPrIhLrfVuyH36dFDNXQoi0+VUyJW997Qu0558zUjK4XM1NOYCW96Uz1ERJQcbJ2aqaqqwpQpU3DllVfiqquuwrJly9Da2iqtoklUh1WalOndkOeMHqM7lWCm+VesK3sSpeW6VgiasaHW2QshIiJb2RpEbrnlFhw7dgxz585FQ0MDLr30Urz66qsRBaxeY7TmIlZLtm4Jm1o4fe4cjpxq0W3+JWdXoHB7Ga9ewJLX4hARUeKzvVj13nvvTbipGL0phHinCATZ3890hhAAEDo/UTmgJOI5sYQg+U1db3ddsWOsXY3bjEqUERsiIrKGp1bNJILHt7wpFVkCsY0g6C251SpKNcvMTV0caXBz3xojxBVJWiuTiIgosTCImKS3QZ2W2a++AgBRC1sB92/4Xh+REOtsuHSXiCg5MIhYRK+uRN4/w8jNfsnWLXjrwP6Ix/WmWGK9TqMjDMrN+eyupyEiouTGIKIQ69SEkboSMx1EDwYClo5MiHUt8Y62iE3XrHo9IiJKbQwiCnZMTczYUKvZQTRWZkci5M3R4hlVKe6VjV3HjmJ/cxNmbKjFqkmTY3odLWa+rse3vIl3Dx1SPZYBiYgoMTCIOEDeJCxe8YxEWNEc7dlJk6U+J1Z+XSIzX9cjY8aZakNPRETewyCSYOQjEYu3bsEmlToSsR18LPUbRp7j1CZ1Wv1WxOsx2nOFiIi8i0EkDkZXv4g3cXGpbyxBQE68GR/SqCMRe4LEMmripSmNtw7s9/QKHiIiih+DSByM1pMob+5Gbvb5WT01a0r0wooYVBZv3YI0mFvqqtwhV66oV7blNSFERJTaGEQUnFqOunjrFgCh5mbym39xr2w8O2myoX4l+Vk9Mb58SETbc3Eprl7jNDXiqphojIzsEBERGcEgouDU1MQhlW3sW9vbsXHf3ojVLcpRCrEGxGxzNau6prrdBl5OrbMtAxIRUeJgEHGJ/GYu/7vaKhC9UQqjAcPrXVNj4aVQRERE5jGIOMhIceuIgkJpJYjR3+ytDBjxFtIaYTQ4OXEtRETkLgYRBxkJDGq/4ce7268ZTkxNxVrkS0REyYdBJA52/cb+mx3b8B+Xj8KSrVuw5+uvw3b7jfe1o7lrQy0EgCtjiIjIMQwicbDrN/Ycvx+A8zvMHpZ1StWbPtHalE9+DBERkR4GEQ/6wfARtryumaXJetMnVm/KR0REqYlBJAkYDRh6oxRFvbIjepIQERHZiUFEh1W9N+xk1XU4VRviVNM4IiLyPgYRHYm2NDYWMzbUSm3hxSZqZgzK7R22+ZxakzE5rwQ4IiJyH4OIg7x6A87KiH06ZlBub7wxZZqFV0NERKmEQcRhyq3t5fKzeppu226FWAJSaU4Ohuf3l0Y/EmEKi4iIvIdBxGFOL8mN1diyQVLIOH3unDR1I04hzRk9JuxrScb28UREZD8GEZcoN7ITFfXKDisaNdIW3g56gUmsK4mlpoSIiEjEIKLDrhUel/Tvj5KcHOnjs+fP4/jp0yjtfEzsrurUSIORqRVxWkk+QkJERBQPBhEddo04TLt8VNTP/4fO563W2n5O9XFB6Pp7okwrERFR4mAQcZDeqIM4XeNGgaeRHiKrd2zDJ0fDR2dKcnLwEwYUIiKKEYOIg/SmWfaedLfgUy8ofXL0KP78+e6wzw3P74+fjB7DJmVERBQTBhGSxFKP0treDsC7PVKIiMjbfG5fACW2/c1NmLGh1u3LICKiBMUREY9TTm3IN6Yrla26cRNX0BARUawYRDzOS1MeJZ3dVOXEqRkiIqJYMIg4yKqCTreanP1k9BjVFTKcmiEiolgxiDhILyAU9coGoB9M3GqnruwGK4aeC/r2dfxaiIgoOTCIeIiRXh520huxEZcXD8rtjayMDKnZGRudERFRrBhESKI3YtPa3o5Bub3xxpRpDl0RERElOwaRFKNVX/LNkhI8MmZcxOPidExrezv2NzdJxapu1akQEVFyYRBJMWbrS7S6vbpVp0JERMmFDc2IiIjINRwRSUDRVtV8s6Qk6nOLemWbGslQnmtIH+4dQ0RE1mEQSUDx1F+smjQZ169bjf3NTbafi4iISA+DSBLS20U3KyPDhasiIiKKZEuNyP79+zFt2jSUl5cjMzMTQ4YMwbx589DOduCOEAtJlf9prXIhIiJyiy0jIrt370YwGMSzzz6LiooK7Ny5E9OnT0draysWL15sxynJBLU6D7Nt5qMdb/a1iIgodaUJgtgf015PPfUUVqxYgS+//NLwcwKBAHJzc9Hc3Iwcj+w0azcr+nNM+t3/qBakDs/vjw3/dkfc10hERBSNmfu3YzUizc3NyNP5TbmtrQ1tbW3Sx4FAwO7L8hz25yAiolTiSB+Ruro6LF++HHfddVfU42pqapCbmyv9V1pa6sTlERERkUtMBZGHHnoIaWlpUf/bvXt32HPq6+vx3e9+FzfffDOmT58e9fWrq6vR3Nws/Xfw4EHzXxERERElDFNTM3PmzMHUqVOjHjN48GDp74cPH8a4ceMwevRorFq1Svf1/X4//H6/mUsiFXq76BIREXmFqSCSn5+P/Px8Q8fW19dj3LhxuOKKK7BmzRr4fOwm7xQ2ISMiokRhS7FqfX09xo4di7KyMixevBjHjh2TPldYWGjHKYmIiCgB2RJENm7ciLq6OtTV1aFEsfeJQ6uFExb7cxARUSpxrI9ILFKxjwgREVGiM3P/ZuEGERERuYZBhIiIiFzDIEJERESuYRAhIiIi1zCIEBERkWsYRIiIiMg1DCJERETkGgYRIiIicg2DCBEREbmGQYSIiIhcY8teM1YRu88HAgGXr4SIiIiMEu/bRnaR8XQQaWlpAQCUlpa6fCVERERkVktLC3Jzc6Me4+lN74LBIA4fPozs7GykpaW5fTm2CwQCKC0txcGDB1N+kz++F134XnThe9GF70UXvhddvPJeCIKAlpYWFBcXw+eLXgXi6RERn8+HkpISty/DcTk5OSn/P5OI70UXvhdd+F504XvRhe9FFy+8F3ojISIWqxIREZFrGESIiIjINQwiHuL3+zFv3jz4/X63L8V1fC+68L3owveiC9+LLnwvuiTie+HpYlUiIiJKbhwRISIiItcwiBAREZFrGESIiIjINQwiRERE5BoGEZcNGjQIaWlpYf898cQTUZ9z9uxZzJo1C3379kWvXr3wL//yL2hsbHToiu2xf/9+TJs2DeXl5cjMzMSQIUMwb948tLe3R33e2LFjI96/u+++26Grts4zzzyDQYMGoUePHqisrMT7778f9fgXX3wRw4YNQ48ePXDJJZfgr3/9q0NXap+amhqMGjUK2dnZ6N+/PyZPnow9e/ZEfc7atWsjvv89evRw6Irt84tf/CLi6xo2bFjU5yTjzwSg/m9kWloaZs2apXp8Mv1MbN68GZMmTUJxcTHS0tJQW1sb9nlBEDB37lwUFRUhMzMT48ePxxdffKH7umb/vbEbg4gHPPbYYzhy5Ij033333Rf1+AceeAAbNmzAiy++iE2bNuHw4cP453/+Z4eu1h67d+9GMBjEs88+i127duHpp5/GypUr8fDDD+s+d/r06WHv36JFixy4Yuv84Q9/QFVVFebNm4cdO3Zg5MiRmDBhAo4ePap6/NatW/Fv//ZvmDZtGj788ENMnjwZkydPxs6dOx2+cmtt2rQJs2bNwrvvvouNGzfi3Llz+M53voPW1taoz8vJyQn7/h84cMChK7bX8OHDw76ut99+W/PYZP2ZAIBt27aFvQ8bN24EANx8882az0mWn4nW1laMHDkSzzzzjOrnFy1ahF/96ldYuXIl3nvvPWRlZWHChAk4e/as5mua/ffGEQK5qqysTHj66acNH9/U1CR0795dePHFF6XHPvvsMwGA8M4779hwhe5ZtGiRUF5eHvWY6667Trj//vuduSCbXHXVVcKsWbOkjzs6OoTi4mKhpqZG9fgf/OAHwo033hj2WGVlpXDXXXfZep1OO3r0qABA2LRpk+Yxa9asEXJzc527KIfMmzdPGDlypOHjU+VnQhAE4f777xeGDBkiBINB1c8n688EAOGll16SPg4Gg0JhYaHw1FNPSY81NTUJfr9f+N3vfqf5Omb/vXECR0Q84IknnkDfvn1x2WWX4amnnsL58+c1j92+fTvOnTuH8ePHS48NGzYMAwcOxDvvvOPE5TqmubkZeXl5use98MIL6NevHy6++GJUV1fj9OnTDlydNdrb27F9+/aw76fP58P48eM1v5/vvPNO2PEAMGHChKT8/gPQ/Rk4deoUysrKUFpaiptuugm7du1y4vJs98UXX6C4uBiDBw/G7bffjq+++krz2FT5mWhvb8fzzz+PO++8M+pGqMn6MyG3b98+NDQ0hH3fc3NzUVlZqfl9j+XfGyd4etO7VPDjH/8Yl19+OfLy8rB161ZUV1fjyJEjWLp0qerxDQ0NyMjIQO/evcMeLygoQENDgwNX7Iy6ujosX74cixcvjnrcbbfdhrKyMhQXF+Pjjz/Ggw8+iD179uBPf/qTQ1can+PHj6OjowMFBQVhjxcUFGD37t2qz2loaFA9Ppm+/8FgELNnz8Y111yDiy++WPO4Cy64AM899xxGjBiB5uZmLF68GKNHj8auXbsSesPMyspKrF27FhdccAGOHDmC+fPnY8yYMdi5cyeys7Mjjk+FnwkAqK2tRVNTE6ZOnap5TLL+TCiJ31sz3/dY/r1xAoOIDR566CE8+eSTUY/57LPPMGzYMFRVVUmPjRgxAhkZGbjrrrtQU1OTUC16tZh5L0T19fX47ne/i5tvvhnTp0+P+twZM2ZIf7/kkktQVFSEG264AXv37sWQIUPiu3hyzaxZs7Bz586odREAcPXVV+Pqq6+WPh49ejQuvPBCPPvss/jP//xPuy/TNhMnTpT+PmLECFRWVqKsrAx//OMfMW3aNBevzF2rV6/GxIkTUVxcrHlMsv5MJDMGERvMmTMnamIHgMGDB6s+XllZifPnz2P//v244IILIj5fWFiI9vZ2NDU1hY2KNDY2orCwMJ7LtoXZ9+Lw4cMYN24cRo8ejVWrVpk+X2VlJYDQiEoiBJF+/fohPT09YtVTtO9nYWGhqeMTzb333ou//OUv2Lx5s+nfYLt3747LLrsMdXV1Nl2dO3r37o2hQ4dqfl3J/jMBAAcOHMD//d//mR7tTNafCfF729jYiKKiIunxxsZGXHrpparPieXfGyewRsQG+fn5GDZsWNT/MjIyVJ/70UcfwefzoX///qqfv+KKK9C9e3e8/vrr0mN79uzBV199FfZbgFeYeS/q6+sxduxYXHHFFVizZg18PvM/nh999BEAhP2P6WUZGRm44oorwr6fwWAQr7/+uub38+qrrw47HgA2btzoye+/GYIg4N5778VLL72EN954A+Xl5aZfo6OjA5988knCfP+NOnXqFPbu3av5dSXrz4TcmjVr0L9/f9x4442mnpesPxPl5eUoLCwM+74HAgG89957mt/3WP69cYRrZbIkbN26VXj66aeFjz76SNi7d6/w/PPPC/n5+cIPf/hD6ZhDhw4JF1xwgfDee+9Jj919993CwIEDhTfeeEP44IMPhKuvvlq4+uqr3fgSLHPo0CGhoqJCuOGGG4RDhw4JR44ckf6THyN/L+rq6oTHHntM+OCDD4R9+/YJL7/8sjB48GDh2muvdevLiMnvf/97we/3C2vXrhU+/fRTYcaMGULv3r2FhoYGQRAE4Y477hAeeugh6fi///3vQrdu3YTFixcLn332mTBv3jyhe/fuwieffOLWl2CJmTNnCrm5ucJbb70V9v0/ffq0dIzyvZg/f77w2muvCXv37hW2b98u3HrrrUKPHj2EXbt2ufElWGbOnDnCW2+9Jezbt0/4+9//LowfP17o16+fcPToUUEQUudnQtTR0SEMHDhQePDBByM+l8w/Ey0tLcKHH34ofPjhhwIAYenSpcKHH34oHDhwQBAEQXjiiSeE3r17Cy+//LLw8ccfCzfddJNQXl4unDlzRnqN66+/Xli+fLn0sd6/N25gEHHR9u3bhcrKSiE3N1fo0aOHcOGFFwoLFy4Uzp49Kx2zb98+AYDw5ptvSo+dOXNGuOeee4Q+ffoIPXv2FP7pn/4p7IadiNasWSMAUP1PpHwvvvrqK+Haa68V8vLyBL/fL1RUVAg//elPhebmZpe+itgtX75cGDhwoJCRkSFcddVVwrvvvit97rrrrhOmTJkSdvwf//hHYejQoUJGRoYwfPhw4ZVXXnH4iq2n9f1fs2aNdIzyvZg9e7b0vhUUFAjf+973hB07djh/8Ra75ZZbhKKiIiEjI0MYMGCAcMsttwh1dXXS51PlZ0L02muvCQCEPXv2RHwumX8m3nzzTdX/J8SvNxgMCj//+c+FgoICwe/3CzfccEPEe1RWVibMmzcv7LFo/964IU0QBMHRIRgiIiKiTqwRISIiItcwiBAREZFrGESIiIjINQwiRERE5BoGESIiInINgwgRERG5hkGEiIiIXMMgQkRERK5hECEiIiLXMIgQERGRaxhEiIiIyDUMIkREROSa/w9hqmm3PtJlSgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# original method, use numpy to implement\n",
    "%matplotlib inline\n",
    "from sklearn.datasets import make_blobs\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import exp, log\n",
    "import copy\n",
    "\n",
    "## Hyper-paramters definition\n",
    "LR = 0.05\n",
    "EPOCH = 20\n",
    "'''\n",
    "numpy.random.normal(loc=0.0, scale=1.0, size=None)  \n",
    "loc:float 概率分布的均值，对应着整个分布的中心center\n",
    "scale:float概率分布的标准差，对应于分布的宽度，scale越大越矮胖，scale越小，越瘦高\n",
    "size:int or tuple of ints\n",
    "我们更经常会用到np.random.randn(size)所谓标准正太分布（μ=0, σ=1），对应于np.random.normal(loc=0, scale=1, size)\n",
    "'''\n",
    "THETA = np.random.normal(0, 0.1, size=(3,3)).reshape(3, 3) # learnable parameters #和上述公式中表述的一样，为列向量组合而成的矩阵\n",
    "# https://www.jianshu.com/p/069d8841bd8e make_blobs函数是为聚类产生数据集\n",
    "X, Y = make_blobs(n_samples=600, centers=3, n_features=2, random_state=3)\n",
    "print(X.shape, Y.shape)\n",
    "# plt.scatter(X[:, 0], X[:, 1], c=Y)\n",
    "# plt.show()\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.1, random_state=0)\n",
    "plt.scatter(X_train[:,0], X_train[:, 1], c=Y_train, edgecolors='white', marker='s')\n",
    "plt.show()\n",
    "\n",
    "X0_train = np.ones([X_train.shape[0],1],dtype=X_train.dtype)\n",
    "X0_test = np.ones([X_test.shape[0],1], dtype=X_test.dtype)\n",
    "X_train_original = copy.deepcopy(X_train)\n",
    "X_train = np.concatenate((X0_train,X_train), axis=1)\n",
    "X_test = np.concatenate((X0_test, X_test), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1/20] loss is: 92.61701043022114\n",
      "[2/20] loss is: 41.38712001322031\n",
      "[3/20] loss is: 32.70813079961782\n",
      "[4/20] loss is: 26.95700278906443\n",
      "[5/20] loss is: 22.719074880715258\n",
      "[6/20] loss is: 19.55606034508767\n",
      "[7/20] loss is: 17.157352074840333\n",
      "[8/20] loss is: 15.304130200953017\n",
      "[9/20] loss is: 13.847470955508832\n",
      "[10/20] loss is: 12.681056335546428\n",
      "[11/20] loss is: 11.727450331954792\n",
      "[12/20] loss is: 10.931356104558194\n",
      "[13/20] loss is: 10.254078410206985\n",
      "[14/20] loss is: 9.668593333017522\n",
      "[15/20] loss is: 9.155759130032143\n",
      "[16/20] loss is: 8.701706332978276\n",
      "[17/20] loss is: 8.296128800377474\n",
      "[18/20] loss is: 7.931177747438058\n",
      "[19/20] loss is: 7.600739370079159\n",
      "[20/20] loss is: 7.299954199099987\n",
      "THETA [[ 3.78130734  0.34763698 -4.0333719 ]\n",
      " [ 0.44964144 -3.13714972  2.78159912]\n",
      " [ 0.62486146 -0.63897881  0.43461416]] (60,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHlhJREFUeJzt3X+MVeWdP/DPzLAM1AwjVEYFBuWHG6jaFEHZSNOVlVQbulvSxm4TTCg16LJYoZgqtFVrW5naZV1a2oC4KZItVrux7ra6aAiNGrM2KFi3VEHxRx0H+eGiMxT7HYQ53z+6TB1hxhmYc5/74/VKbuKcc6/nczJm7tvn+TzPqcqyLAsAgASqUxcAAFQuQQQASEYQAQCSEUQAgGQEEQAgGUEEAEhGEAEAkhFEAIBkBqQuoCcdHR2xa9euqKuri6qqqtTlAAC9kGVZHDhwIEaMGBHV1T2PeRR1ENm1a1c0NjamLgMAOAHNzc0xatSoHt9T1EGkrq4uIv50I0OGDElcDQDQG21tbdHY2Nj5Pd6Tog4iR6djhgwZIogAQInpTVuFZlUAIBlBBABIRhABAJIRRACAZAQRACAZQQQASEYQAQCSEUQAgGQEEQAgGUEEAEimqLd4BwD6V8fb10ccfun4JweMi+pT/7mg9QgiAFAmegwZtZ+I6rrFfzp/+LnCFtYDQQSAitHtF3WCkYBc9BQyasYWtpZeEkQAqBxFNhqAZlUAICFBBABIRhABAJLRIwIA5WLAuB7OjerFe3o4lxNBBADKRG9W/hTb6iBBBIDK0d3/8ScYCeBPBBEAKkaxjQagWRUASEgQAQCSEUQAgGQEEQAgGUEEAEhGEAEAkhFEAIBkBBEAIBlBBABIxs6qAJBYx9vXRxx+6fgnB4wr6x1hBREASO3wSxGHn0tdRRKmZgCAZAQRACAZQQQASEYQAQCS0awKAKkNGHdi5z5AKazGEUQAILHcAkEJrMYxNQMAJCOIAADJCCIAQDKCCACQjGZVAChXOa3G6U+CCACUqWJYnvtBTM0AAMkYEQGAItbtpmRFsiHZyRJEACgppbBbaL8qgU3JToYgAkBpKfMv5kqjRwQASEYQAQCSyS2IHDlyJG666aYYM2ZMDB48OMaNGxff/va3I8uyvC4JQCWrOSuiZmzqKuij3HpEbr/99li1alWsW7cuzj333Hj66adj7ty5UV9fH9ddd11elwWgEtWcFdXDN6auIh/dbTxWJBuSnazcgsh///d/x2c+85mYOXNmREScffbZ8dOf/jQ2b96c1yUBqFRVp0REeS51LdW6eyu3IHLxxRfHmjVr4oUXXoi//Mu/jGeffTaeeOKJuOOOO7r9THt7e7S3t3f+3NbWlld5AJSq440EHJ2SsaKm5OQWRJYsWRJtbW0xYcKEqKmpiSNHjsRtt90Ws2fP7vYzTU1Nceutt+ZVEgBloNxHCCpNbs2qP/vZz2L9+vVxzz33xNatW2PdunWxfPnyWLduXbefWbp0abS2tna+mpub8yoPACgCuY2IfPWrX40lS5bEF77whYiIOP/88+P3v/99NDU1xZw5c477mdra2qitrc2rJACgyOQWRN55552oru464FJTUxMdHR15XRIAKkI5NeXmFkT+9m//Nm677bYYPXp0nHvuufHMM8/EHXfcEV/60pfyuiQAla7Ml7p26kVTbqmEldyCyMqVK+Omm26Kf/zHf4y9e/fGiBEj4pprrombb745r0sCUOGK6Qu2P3WGipqxUT20+9WnXZTICqLcgkhdXV2sWLEiVqxYkdclAKAylEioOBGeNQMAJCOIAADJ5DY1AwDkpIyacgURACgV2cGI6GVTbomEFUEEAIrde8JDx/75EVWD33NuVFTXLT7mI6WygkgQAYAiVyqh4kRoVgUAkhFEAIBkBBEAIBlBBABIRhABAJIRRACAZAQRACAZQQQASMaGZgDwfzrevj7i8EvHnhgwrqw3FUtJEAGAow6/FHH4udRVVBRTMwBAMoIIAJCMIAIAJCOIAADJaFYFgKMGjOvbcU6aIAIA/8cS3cIzNQMAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDIDUhcAUGyarvx+vPZ8y3HPjZ44Mpb+ZGGBK4LyJYgAvM9rz7fEzmdeSV0GVARTMwBAMoIIAJCMIAIAJKNHBIiI7hs0NWcCeRJEgIjQoPleoyeOPKFzQN8JIkASxTwCk/r6UEkEESAJIzBAhGZVACAhIyLAMUaOPyMG1w2OiIjGCXoigPwIIkBE/LkJc3jjafGt/7ghcTVApRBEoAL0pjH0/Q2axdxMCpSPXINIS0tL3HjjjbFhw4Z45513Yvz48bF27dqYMmVKnpcF3udEGkM1kwKFkFsQeeutt2LatGkxffr02LBhQwwfPjxefPHFGDp0aF6XBEpId/tx2KeDctbx9vURh1869sSAcVF96j8XvqAikFsQuf3226OxsTHWrl3beWzMmDF5XQ4oMaZ3qEiHX4o4/FzqKopKbst3f/GLX8SUKVPiiiuuiIaGhpg0aVLcddddeV0OAChBuQWRl19+OVatWhXnnHNOPPLIIzF//vy47rrrYt26dd1+pr29Pdra2rq8AIDyldvUTEdHR0yZMiWWLVsWERGTJk2Kbdu2xerVq2POnDnH/UxTU1PceuuteZUEFetE+jH0cACFkFsQOfPMM+MjH/lIl2MTJ06M+++/v9vPLF26NBYvXtz5c1tbWzQ2NuZVIlSME+nH0MMBFEJuQWTatGmxY8eOLsdeeOGFOOuss7r9TG1tbdTW1uZVEgCkNWBc345XgNyCyFe+8pW4+OKLY9myZfH5z38+Nm/eHGvWrIk1a9bkdUngfWxKBsWlUpfo9iS3IHLhhRfGAw88EEuXLo1vfetbMWbMmFixYkXMnj07r0sC72NTMqDY5bqz6qc//en49Kc/neclAIASltvyXQCADyKIAADJCCIAQDK59ogAadmUDCh2ggiUMUt0gWJnagYASMaICNArNkcD8iCIAL1iczQgD6ZmAIBkBBEAIBlBBABIRhABAJIRRACAZKyaAXrFLq1AHgQRoFfsFQLkwdQMAJCMIAIAJCOIAADJCCIAQDKaVQEPtAOSEUQAD7QDkjE1AwAkI4gAAMkIIgBAMoIIAJCMZlXAc2SAZAQRwBJdIBlTMwBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwt3oGkmq78frz2fMsxx0dPHGnreagAggiQ1GvPt8TOZ15JXQaQiCACJcYIAlBOBBEoMUYQgHKiWRUASMaICNDtdE+EKR8gX4IIkHS6Z/TEkX06DpQXQQRIymgLVDZBBIrQ8aZKGieMjK+tX2gEASgrgggUoeNNlfzxwB8jwggCUF6smoES0bJzd9w863upywDoV0ZEoITsa36zV+/r66ZnPU3rmPIB8iSIQBnq6yoY0z1AKgWbmvnud78bVVVVsWjRokJdEgAocgUZEXnqqafizjvvjI9+9KOFuByUvJNdGdM44fjvG9542gnXBJCH3IPIH/7wh5g9e3bcdddd8Z3vfCfvy0FZONmpkq+tN9UClIbcg8iCBQti5syZMWPGDEEEcmardqDU5BpE7r333ti6dWs89dRTvXp/e3t7tLe3d/7c1taWV2lQNpqu/H5k2Z9GQYrxybx9XcEDVJbcgkhzc3MsXLgwNm7cGIMGDerVZ5qamuLWW2/NqyQoS92NgBSLYgxHQPHIbdXMli1bYu/evXHBBRfEgAEDYsCAAfHYY4/FD37wgxgwYEAcOXLkmM8sXbo0WltbO1/Nzc15lQcAFIHcRkQuvfTS+O1vf9vl2Ny5c2PChAlx4403Rk1NzTGfqa2tjdra2rxKAgCKTG5BpK6uLs4777wux0455ZT48Ic/fMxxAKAy2VkV6LOTaUAdOf6MGFw3uNu9ToDKUtAg8uijjxbyclARRk8cGVn253/u6X39pS8NqO+97vDG0+Jb/3FDv9UBlD4jIlDi3jsCkedy2PcuE+6L49VkSS9wlCAC9Ep/LhO2pBc4qmAPvQMAeD9BBABIxtQM0Gcn+3RggKMEEaDX/njgjxGRb1MsUFkEEaBXRk8cGa893xKrFq+L+XfMiWWzvx/N249tYO3NEl0jKsBRggjQK+8fBWnefuIrX4yoAEcJIlDG7NcBFDtBBMqY/TqAYmf5LgCQjBER4IRoOAX6gyACnBA9JkB/EESA3GmaBbojiEAZK5bpE02zQHcEEShjRhuAYmfVDACQjCACACRjagbKnEZRoJgJIlDmiqFRtFiaZoHiI4gAuTPyAnRHjwgAkIwgAgAkI4gAAMnoEYEyp1EUKGaCCJQ5jaJAMTM1AwAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQjCACACRjZ1WgW01Xfj9ee77luOdGTxxp11bgpAkiQLdee74ldj7zSuoygDJmagYASEYQAQCSMTVTgr7yyEOxc//+Y46PHzYs/uWymQkqoph019ehpwMoRoJICdq5f3/8bt/e1GVQpFL3dWhwBfpCEAG6NXriyD6fSx2EgNIiiADdMnoB5E2zKgCQjCACACRjaqYEjR82rE/HqSzd9W701O8BkIogUoIs0aUn+jqAUiKIQAX48Tfuid2v7Iv/d7A99jW/2eVcfy+pPZGVNkDlEkSgAjy14TcFW1JrRAboC82qAEAygggAkEyuQaSpqSkuvPDCqKuri4aGhpg1a1bs2LEjz0sCACUk1yDy2GOPxYIFC+LXv/51bNy4Md5999345Cc/GQcPHszzsgBAici1WfXhhx/u8vPdd98dDQ0NsWXLlvjEJz6R56UBgBJQ0FUzra2tERExrJuNt9rb26O9vb3z57a2toLUBeXOklqgWFVlWZYV4kIdHR3xd3/3d/H222/HE088cdz3fPOb34xbb731mOOtra0xZMiQvEsEAPpBW1tb1NfX9+r7u2BBZP78+bFhw4Z44oknYtSoUcd9z/FGRBobGwURACghfQkiBZmaufbaa+PBBx+Mxx9/vNsQEhFRW1sbtbW1hSip7H3lkYdi5/79xxwfP2yYLeIBKBq5BpEsy+LLX/5yPPDAA/Hoo4/GmDFj8rwc77Fz//743b69qcuAktF05ffjtedbjjne31vgA13lGkQWLFgQ99xzT/znf/5n1NXVxe7duyMior6+PgYPHpznpQH65LXnWwq2DT7wZ7nuI7Jq1apobW2NSy65JM4888zO13333ZfnZQGAEpH71AwAQHc8awYASEYQAQCSKejOqhTO+G52r+3uOFS67naYtfMs5KtgG5qdiL5siAIAFIe+fH+bmgEAkhFEAIBkBBEAIBnNqkWou+fERHhWDADlRRApQp4TA0ClMDUDACQjiAAAyQgiAEAygggAkIxm1SLU0zbsR88teviheOmt/VbRAFDSBJEi1Jtg8T97dserrW/nXwwA5EgQOQmF3O/j6AhIRMTBQ4eEEADKgiByEgq538dLb9lbBIDyo1kVAEhGEAEAkjE1UyK6W0nT0wobACh2gkiJsEQXgHIkiJyE3uz3AQB0TxA5CUYpAODkaFYFAJIRRACAZAQRACAZQQQASEYQAQCSEUQAgGQs3y0D3T0FuL+fAAwA/U0QKQOFfAowAPQnUzMAQDKCCACQjCACACQjiAAAyWhWLQPdPenXE4ABKHaCSBkolyW63S1DjrAUGaBcCSIUDcuQASqPHhEAIBkjIkXMVAUA5U4QKWKmKgAod6ZmAIBkBBEAIBlTMxSNnvY9sScKQHkSRCgamm8BKo8gUsQKOUJghQ4AKQgiRayQX/5W6ACQgmZVACAZQQQASEYQAQCSyT2I/OhHP4qzzz47Bg0aFFOnTo3NmzfnfUkAoETk2qx63333xeLFi2P16tUxderUWLFiRVx22WWxY8eOaGhoyPPSuepuhUmK1SX9VYs9PABIIdcgcscdd8S8efNi7ty5ERGxevXqeOihh+LHP/5xLFmyJM9L56qYVpj0Vy2W5wKQQm5TM4cOHYotW7bEjBkz/nyx6uqYMWNGPPnkk8f9THt7e7S1tXV5AQDlK7cg8uabb8aRI0fi9NNP73L89NNPj927dx/3M01NTVFfX9/5amxszKs8AKAIFNWqmaVLl0Zra2vnq7m5OXVJAECOcusROe2006Kmpib27NnT5fiePXvijDPOOO5namtro7a2Nq+SAIAik1sQGThwYEyePDk2bdoUs2bNioiIjo6O2LRpU1x77bV5XbbidLeixUoXAEpBrqtmFi9eHHPmzIkpU6bERRddFCtWrIiDBw92rqIpVcX05W+1CwClLNcg8vd///exb9++uPnmm2P37t3xsY99LB5++OFjGlhLjS9/AOgfVVmWZamL6E5bW1vU19dHa2trDBkyJHU5AEAv9OX7u6hWzQAAlUUQAQCSEUQAgGRybVblg139y/+IN/5woPPnE31wXjE9iA8AeksQydkHBYQ3/nCgXx5aV0wP4gOA3hJEciYgAED39IgAAMkIIgBAMoIIAJCMHpHExg3t+nyaE31eTTE9/wYAeksQydkHBYQVl/fP0lpLdAEoRRUVRLpbShuR334bAgIAdK+igoiltABQXCoqiHDi7NwKQB4EEXrlvaNJZ9efGqcMHBgREVmWsioASp0gQp+cXX9q/GrOVanLAKBMCCL0ydGRkBSNvwCUn4oKIj3tqWG/jb7R+AtAf6ioIOL/0gGguNjiHQBIpqJGRDhxR6eu3r8lPQCcDEGEXjGtBUAeBBFOiMZfAPqDIMIJMUICQH/QrAoAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQzIDUBVSSrzzyUOzcv/+458YPGxb/ctnMAlcEAGkJIgW0c//++N2+vanLAICiYWoGAEhGEAEAkhFEAIBkBBEAIBlBBABIxqqZAho/bNgJnQOAciWIFJB9QgCgK1MzAEAyuQSRV199Na666qoYM2ZMDB48OMaNGxe33HJLHDp0KI/LAQAlKpepme3bt0dHR0fceeedMX78+Ni2bVvMmzcvDh48GMuXL8/jkgBACarKsiwrxIX+6Z/+KVatWhUvv/xyrz/T1tYW9fX10draGkOGDMmxOgCgv/Tl+7tgzaqtra0x7ANWhrS3t0d7e3vnz21tbXmXBQAkVJBm1Z07d8bKlSvjmmuu6fF9TU1NUV9f3/lqbGwsRHkAQCJ9CiJLliyJqqqqHl/bt2/v8pmWlpa4/PLL44orroh58+b1+O9funRptLa2dr6am5v7fkcAQMnoU4/Ivn374n//9397fM/YsWNj4MCBERGxa9euuOSSS+Kv/uqv4u67747q6r4NwOgRAYDSk1uPyPDhw2P48OG9em9LS0tMnz49Jk+eHGvXru1zCAEAyl8uzaotLS1xySWXxFlnnRXLly+Pffv2dZ4744wz8rgkAFCCcgkiGzdujJ07d8bOnTtj1KhRXc71ZbXw0fdaPQMApePo93ZvvvMLto/IiXj99detnAGAEtXc3HzMgMT7FXUQ6ejoiF27dkVdXV1UVVV1+762trZobGyM5ubmimtqreR7j6js+3fv7t29V45Su/csy+LAgQMxYsSID+wRLeqn71ZXV39gknqvIUOGlMQvKA+VfO8RlX3/7t29Vxr3Xhr3Xl9f36v3WcoCACQjiAAAyZRFEKmtrY1bbrklamtrU5dScJV87xGVff/u3b1XGvdenvde1M2qAEB5K4sREQCgNAkiAEAygggAkIwgAgAkU7ZB5KGHHoqpU6fG4MGDY+jQoTFr1qzUJRVUe3t7fOxjH4uqqqr4zW9+k7qc3L366qtx1VVXxZgxY2Lw4MExbty4uOWWW+LQoUOpS8vFj370ozj77LNj0KBBMXXq1Ni8eXPqknLX1NQUF154YdTV1UVDQ0PMmjUrduzYkbqsJL773e9GVVVVLFq0KHUpBdHS0hJXXnllfPjDH47BgwfH+eefH08//XTqsnJ35MiRuOmmm7r8Xfv2t7/dp2e2lYKi3ln1RN1///0xb968WLZsWfzN3/xNHD58OLZt25a6rIK64YYbYsSIEfHss8+mLqUgtm/fHh0dHXHnnXfG+PHjY9u2bTFv3rw4ePBgLF++PHV5/eq+++6LxYsXx+rVq2Pq1KmxYsWKuOyyy2LHjh3R0NCQurzcPPbYY7FgwYK48MIL4/Dhw/G1r30tPvnJT8Zzzz0Xp5xySuryCuapp56KO++8Mz760Y+mLqUg3nrrrZg2bVpMnz49NmzYEMOHD48XX3wxhg4dmrq03N1+++2xatWqWLduXZx77rnx9NNPx9y5c6O+vj6uu+661OX1n6zMvPvuu9nIkSOzf/3Xf01dSjL/9V//lU2YMCH73e9+l0VE9swzz6QuKYnvfe972ZgxY1KX0e8uuuiibMGCBZ0/HzlyJBsxYkTW1NSUsKrC27t3bxYR2WOPPZa6lII5cOBAds4552QbN27M/vqv/zpbuHBh6pJyd+ONN2Yf//jHU5eRxMyZM7MvfelLXY599rOfzWbPnp2oonyU3dTM1q1bo6WlJaqrq2PSpElx5plnxqc+9amKGRHZs2dPzJs3L/7t3/4tPvShD6UuJ6nW1tYYNmxY6jL61aFDh2LLli0xY8aMzmPV1dUxY8aMePLJJxNWVnitra0REWX3O+7JggULYubMmV1+/+XuF7/4RUyZMiWuuOKKaGhoiEmTJsVdd92VuqyCuPjii2PTpk3xwgsvRETEs88+G0888UR86lOfSlxZ/yq7IPLyyy9HRMQ3v/nN+MY3vhEPPvhgDB06NC655JLYv39/4urylWVZfPGLX4x/+Id/iClTpqQuJ6mdO3fGypUr45prrkldSr96880348iRI3H66ad3OX766afH7t27E1VVeB0dHbFo0aKYNm1anHfeeanLKYh77703tm7dGk1NTalLKaiXX345Vq1aFeecc0488sgjMX/+/Ljuuuti3bp1qUvL3ZIlS+ILX/hCTJgwIf7iL/4iJk2aFIsWLYrZs2enLq1flUwQWbJkSVRVVfX4OtonEBHx9a9/PT73uc/F5MmTY+3atVFVVRX//u//nvguTkxv733lypVx4MCBWLp0aeqS+01v7/29Wlpa4vLLL48rrrgi5s2bl6hy8rRgwYLYtm1b3HvvvalLKYjm5uZYuHBhrF+/PgYNGpS6nILq6OiICy64IJYtWxaTJk2Kq6++OubNmxerV69OXVrufvazn8X69evjnnvuia1bt8a6deti+fLlZRfCSqZZ9frrr48vfvGLPb5n7Nix8cYbb0RExEc+8pHO47W1tTF27Nh47bXX8iwxN72991/96lfx5JNPHvMsgilTpsTs2bNL8j/e3t77Ubt27Yrp06fHxRdfHGvWrMm5usI77bTToqamJvbs2dPl+J49e+KMM85IVFVhXXvttfHggw/G448/HqNGjUpdTkFs2bIl9u7dGxdccEHnsSNHjsTjjz8eP/zhD6O9vT1qamoSVpifM888s8vf84iIiRMnxv3335+oosL56le/2jkqEhFx/vnnx+9///toamqKOXPmJK6u/5RMEBk+fHgMHz78A983efLkqK2tjR07dsTHP/7xiIh4991349VXX42zzjor7zJz0dt7/8EPfhDf+c53On/etWtXXHbZZXHffffF1KlT8ywxN72994g/jYRMnz69cxSsurpkBvx6beDAgTF58uTYtGlT55L0jo6O2LRpU1x77bVpi8tZlmXx5S9/OR544IF49NFHY8yYMalLKphLL700fvvb33Y5Nnfu3JgwYULceOONZRtCIiKmTZt2zDLtF154oWT/nvfFO++8c8zfsZqams6R/7KRuls2DwsXLsxGjhyZPfLII9n27duzq666KmtoaMj279+furSCeuWVVypm1czrr7+ejR8/Prv00kuz119/PXvjjTc6X+Xm3nvvzWpra7O77747e+6557Krr746O/XUU7Pdu3enLi1X8+fPz+rr67NHH320y+/3nXfeSV1aEpWyambz5s3ZgAEDsttuuy178cUXs/Xr12cf+tCHsp/85CepS8vdnDlzspEjR2YPPvhg9sorr2Q///nPs9NOOy274YYbUpfWr8oyiBw6dCi7/vrrs4aGhqyuri6bMWNGtm3bttRlFVwlBZG1a9dmEXHcVzlauXJlNnr06GzgwIHZRRddlP36179OXVLuuvv9rl27NnVpSVRKEMmyLPvlL3+ZnXfeeVltbW02YcKEbM2aNalLKoi2trZs4cKF2ejRo7NBgwZlY8eOzb7+9a9n7e3tqUvrV1VZVmZbtAEAJaP8JtEBgJIhiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQzP8H56ntsE6Qs5IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## 2. Make a hypothesis\n",
    "# theta (3,3) x(3，1)\n",
    "import numpy as np\n",
    "\n",
    "# Softmax概率计算，基于SGD单个样本的计算\n",
    "def hypothesis(x, THETA):\n",
    "    x = np.reshape(x, (3, 1))\n",
    "    temp = np.matmul(THETA.T, x)\n",
    "    temp = np.exp(temp)\n",
    "    denominator = np.sum(temp)\n",
    "    hypothesis = temp / denominator # normalize into 1 return \n",
    "    return hypothesis\n",
    "\n",
    "## 3. Loss definition，h_x的shape为[3,1]\n",
    "def compute_loss(x, y,THETA):\n",
    "    loss = 0\n",
    "    x = np.reshape(x, (3, 1))\n",
    "    h_x = hypothesis(x, THETA)    #hypothesis (3, 1)\n",
    "    label = y\n",
    "    loss += (-np.log(h_x[label][0] + 0.0000001))  # loss = - 1{y=k} * log(p(y))\n",
    "    #print(f'loss = {loss}')\n",
    "    return loss\n",
    "\n",
    "## 4. Parameters updating\n",
    "def update_parameters(THETA, x, y):\n",
    "    x = np.reshape(x, (3, 1))\n",
    "    #y = np.reshape(y, (3, 1))\n",
    "    \n",
    "    h_x = hypothesis(x, THETA)\n",
    "    label = y\n",
    "    # 参数更新规则具体详见ppt中的说明或上面的公式，由于梯度有负号，更新规则负负得正，所以在此处更新时表达式后面为正号\n",
    "    for k in range(3):\n",
    "        indicator_val = 1 if label == k else 0\n",
    "        THETA[:, k] = THETA[:, k]+LR*(indicator_val-h_x[k][0])*x[:,0]\n",
    "    return THETA\n",
    "\n",
    "for epoch in range(EPOCH):\n",
    "    #LR = LR * (1 / (1 + DECAY_RATE * epoch))\n",
    "    i = 0 # retrieve H_x\n",
    "    loss = 0\n",
    "    for x,y in zip(X_train,Y_train):\n",
    "        loss += compute_loss(x, y, THETA)\n",
    "        #print('[{0}/{1}] loss is: {2}'.format(epoch+1, EPOCH, loss))\n",
    "        THETA = update_parameters(THETA, x, y)\n",
    "    print('[{0}/{1}] loss is: {2}'.format(epoch+1, EPOCH, loss))\n",
    "\n",
    "i = 0\n",
    "print('THETA', THETA, Y_test.shape)\n",
    "H_test = np.zeros((X_test.shape[0], X_test.shape[1]))\n",
    "#H_test = np.zeros([Y_test.shape[0], 1], dtype=Y_test.dtype)\n",
    "for x, y in zip(X_test, Y_test):\n",
    "    H_test[i] = hypothesis(x, THETA).T\n",
    "    i+=1\n",
    "plt.figure(1)\n",
    "x = np.linspace(-7, 4, 50)\n",
    "plt.scatter(X_test[:, 1], X_test[:, 2], c=np.argmax(H_test, axis=1), edgecolors='white', marker='s')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(540, 3)\n"
     ]
    }
   ],
   "source": [
    "print(X_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1/100] loss is: 158.64259649453635\n",
      "[2/100] loss is: 43.79754199088418\n",
      "[3/100] loss is: 86.37980188219522\n",
      "[4/100] loss is: 46.855269150454866\n",
      "[5/100] loss is: 29.922506324699057\n",
      "[6/100] loss is: 43.476901505766065\n",
      "[7/100] loss is: 23.38222945552819\n",
      "[8/100] loss is: 28.953552369448623\n",
      "[9/100] loss is: 10.240980875473587\n",
      "[10/100] loss is: 7.706053004773441\n",
      "[11/100] loss is: 11.11730251114659\n",
      "[12/100] loss is: 8.009225666764406\n",
      "[13/100] loss is: 8.417449733554589\n",
      "[14/100] loss is: 15.87827815137313\n",
      "[15/100] loss is: 6.223466858744404\n",
      "[16/100] loss is: 9.363186211360794\n",
      "[17/100] loss is: 5.709317741078912\n",
      "[18/100] loss is: 5.353407985550361\n",
      "[19/100] loss is: 6.71406004704416\n",
      "[20/100] loss is: 5.404411690035232\n",
      "[21/100] loss is: 6.1063274050152625\n",
      "[22/100] loss is: 5.267000544169731\n",
      "[23/100] loss is: 4.690486533054402\n",
      "[24/100] loss is: 4.210001575764718\n",
      "[25/100] loss is: 4.922837012948339\n",
      "[26/100] loss is: 4.346951350829265\n",
      "[27/100] loss is: 4.756713860673975\n",
      "[28/100] loss is: 3.5896814730760407\n",
      "[29/100] loss is: 3.4026472192940354\n",
      "[30/100] loss is: 5.015917416792035\n",
      "[31/100] loss is: 3.85889638384758\n",
      "[32/100] loss is: 3.088163253868977\n",
      "[33/100] loss is: 2.936968646940425\n",
      "[34/100] loss is: 4.0552554114919275\n",
      "[35/100] loss is: 5.12947506241476\n",
      "[36/100] loss is: 3.41343946054755\n",
      "[37/100] loss is: 3.2947213457914994\n",
      "[38/100] loss is: 3.164059683236951\n",
      "[39/100] loss is: 2.878916600823181\n",
      "[40/100] loss is: 2.995240223926511\n",
      "[41/100] loss is: 2.3233473718279782\n",
      "[42/100] loss is: 2.3987209070663496\n",
      "[43/100] loss is: 2.97278618665698\n",
      "[44/100] loss is: 2.940057328959757\n",
      "[45/100] loss is: 2.2377564491251736\n",
      "[46/100] loss is: 2.417414661607187\n",
      "[47/100] loss is: 2.6033666454201803\n",
      "[48/100] loss is: 2.535638499096362\n",
      "[49/100] loss is: 2.8563089357361764\n",
      "[50/100] loss is: 2.6023267478852885\n",
      "[51/100] loss is: 2.5359717703381257\n",
      "[52/100] loss is: 3.0105325084602432\n",
      "[53/100] loss is: 2.363325764500458\n",
      "[54/100] loss is: 2.6351886775002327\n",
      "[55/100] loss is: 2.06840882011064\n",
      "[56/100] loss is: 2.0615414095962223\n",
      "[57/100] loss is: 2.186673954090291\n",
      "[58/100] loss is: 2.3320889771730884\n",
      "[59/100] loss is: 2.0164388047019663\n",
      "[60/100] loss is: 2.537690265183396\n",
      "[61/100] loss is: 1.9584891782952176\n",
      "[62/100] loss is: 2.3796866153105976\n",
      "[63/100] loss is: 2.1883566771028042\n",
      "[64/100] loss is: 2.1906916032817407\n",
      "[65/100] loss is: 2.438594445652849\n",
      "[66/100] loss is: 1.9396693217051015\n",
      "[67/100] loss is: 2.2999863481282627\n",
      "[68/100] loss is: 2.0294077727602535\n",
      "[69/100] loss is: 2.343990919979844\n",
      "[70/100] loss is: 2.2790810331204767\n",
      "[71/100] loss is: 2.3345366667186878\n",
      "[72/100] loss is: 2.492956523774057\n",
      "[73/100] loss is: 2.0677993738723\n",
      "[74/100] loss is: 1.981539828003042\n",
      "[75/100] loss is: 1.805640844222927\n",
      "[76/100] loss is: 1.797263601611596\n",
      "[77/100] loss is: 1.9657484772997817\n",
      "[78/100] loss is: 1.8543730247875023\n",
      "[79/100] loss is: 1.9292102934208961\n",
      "[80/100] loss is: 1.839248509724119\n",
      "[81/100] loss is: 2.0508270413205403\n",
      "[82/100] loss is: 1.7589316373071382\n",
      "[83/100] loss is: 1.763465329009667\n",
      "[84/100] loss is: 2.19418236444198\n",
      "[85/100] loss is: 2.127237986691216\n",
      "[86/100] loss is: 2.051355343423311\n",
      "[87/100] loss is: 1.6477540445484793\n",
      "[88/100] loss is: 1.726232601294217\n",
      "[89/100] loss is: 2.078325081439886\n",
      "[90/100] loss is: 1.921286620087866\n",
      "[91/100] loss is: 1.5651126991487507\n",
      "[92/100] loss is: 1.8361334788846482\n",
      "[93/100] loss is: 1.7021831110531744\n",
      "[94/100] loss is: 1.8586933815655844\n",
      "[95/100] loss is: 1.6445705046579595\n",
      "[96/100] loss is: 1.8969868916671724\n",
      "[97/100] loss is: 1.9146952576662941\n",
      "[98/100] loss is: 1.5085549000665748\n",
      "[99/100] loss is: 1.7185195344396778\n",
      "[100/100] loss is: 1.640246021307607\n",
      "THETA [[ 6.6452891   0.57547759 -7.26369884]\n",
      " [ 1.04566343 -4.72919835  3.70367629]\n",
      " [ 0.50815069 -1.52501029  1.1217431 ]] (60,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHlhJREFUeJzt3X+MVeWdP/DPzLAM1AwjVEYFBuWHG6jaFEHZSNOVlVQbulvSxm4TTCg16LJYoZgqtFVrW5naZV1a2oC4KZItVrux7ra6aAiNGrM2KFi3VEHxRx0H+eGiMxT7HYQ53z+6TB1hxhmYc5/74/VKbuKcc6/nczJm7tvn+TzPqcqyLAsAgASqUxcAAFQuQQQASEYQAQCSEUQAgGQEEQAgGUEEAEhGEAEAkhFEAIBkBqQuoCcdHR2xa9euqKuri6qqqtTlAAC9kGVZHDhwIEaMGBHV1T2PeRR1ENm1a1c0NjamLgMAOAHNzc0xatSoHt9T1EGkrq4uIv50I0OGDElcDQDQG21tbdHY2Nj5Pd6Tog4iR6djhgwZIogAQInpTVuFZlUAIBlBBABIRhABAJIRRACAZAQRACAZQQQASEYQAQCSEUQAgGQEEQAgGUEEAEimqLd4BwD6V8fb10ccfun4JweMi+pT/7mg9QgiAFAmegwZtZ+I6rrFfzp/+LnCFtYDQQSAitHtF3WCkYBc9BQyasYWtpZeEkQAqBxFNhqAZlUAICFBBABIRhABAJLRIwIA5WLAuB7OjerFe3o4lxNBBADKRG9W/hTb6iBBBIDK0d3/8ScYCeBPBBEAKkaxjQagWRUASEgQAQCSEUQAgGQEEQAgGUEEAEhGEAEAkhFEAIBkBBEAIBlBBABIxs6qAJBYx9vXRxx+6fgnB4wr6x1hBREASO3wSxGHn0tdRRKmZgCAZAQRACAZQQQASEYQAQCS0awKAKkNGHdi5z5AKazGEUQAILHcAkEJrMYxNQMAJCOIAADJCCIAQDKCCACQjGZVAChXOa3G6U+CCACUqWJYnvtBTM0AAMkYEQGAItbtpmRFsiHZyRJEACgppbBbaL8qgU3JToYgAkBpKfMv5kqjRwQASEYQAQCSyS2IHDlyJG666aYYM2ZMDB48OMaNGxff/va3I8uyvC4JQCWrOSuiZmzqKuij3HpEbr/99li1alWsW7cuzj333Hj66adj7ty5UV9fH9ddd11elwWgEtWcFdXDN6auIh/dbTxWJBuSnazcgsh///d/x2c+85mYOXNmREScffbZ8dOf/jQ2b96c1yUBqFRVp0REeS51LdW6eyu3IHLxxRfHmjVr4oUXXoi//Mu/jGeffTaeeOKJuOOOO7r9THt7e7S3t3f+3NbWlld5AJSq440EHJ2SsaKm5OQWRJYsWRJtbW0xYcKEqKmpiSNHjsRtt90Ws2fP7vYzTU1Nceutt+ZVEgBloNxHCCpNbs2qP/vZz2L9+vVxzz33xNatW2PdunWxfPnyWLduXbefWbp0abS2tna+mpub8yoPACgCuY2IfPWrX40lS5bEF77whYiIOP/88+P3v/99NDU1xZw5c477mdra2qitrc2rJACgyOQWRN55552oru464FJTUxMdHR15XRIAKkI5NeXmFkT+9m//Nm677bYYPXp0nHvuufHMM8/EHXfcEV/60pfyuiQAla7Ml7p26kVTbqmEldyCyMqVK+Omm26Kf/zHf4y9e/fGiBEj4pprrombb745r0sCUOGK6Qu2P3WGipqxUT20+9WnXZTICqLcgkhdXV2sWLEiVqxYkdclAKAylEioOBGeNQMAJCOIAADJ5DY1AwDkpIyacgURACgV2cGI6GVTbomEFUEEAIrde8JDx/75EVWD33NuVFTXLT7mI6WygkgQAYAiVyqh4kRoVgUAkhFEAIBkBBEAIBlBBABIRhABAJIRRACAZAQRACAZQQQASMaGZgDwfzrevj7i8EvHnhgwrqw3FUtJEAGAow6/FHH4udRVVBRTMwBAMoIIAJCMIAIAJCOIAADJaFYFgKMGjOvbcU6aIAIA/8cS3cIzNQMAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDIDUhcAUGyarvx+vPZ8y3HPjZ44Mpb+ZGGBK4LyJYgAvM9rz7fEzmdeSV0GVARTMwBAMoIIAJCMIAIAJKNHBIiI7hs0NWcCeRJEgIjQoPleoyeOPKFzQN8JIkASxTwCk/r6UEkEESAJIzBAhGZVACAhIyLAMUaOPyMG1w2OiIjGCXoigPwIIkBE/LkJc3jjafGt/7ghcTVApRBEoAL0pjH0/Q2axdxMCpSPXINIS0tL3HjjjbFhw4Z45513Yvz48bF27dqYMmVKnpcF3udEGkM1kwKFkFsQeeutt2LatGkxffr02LBhQwwfPjxefPHFGDp0aF6XBEpId/tx2KeDctbx9vURh1869sSAcVF96j8XvqAikFsQuf3226OxsTHWrl3beWzMmDF5XQ4oMaZ3qEiHX4o4/FzqKopKbst3f/GLX8SUKVPiiiuuiIaGhpg0aVLcddddeV0OAChBuQWRl19+OVatWhXnnHNOPPLIIzF//vy47rrrYt26dd1+pr29Pdra2rq8AIDyldvUTEdHR0yZMiWWLVsWERGTJk2Kbdu2xerVq2POnDnH/UxTU1PceuuteZUEFetE+jH0cACFkFsQOfPMM+MjH/lIl2MTJ06M+++/v9vPLF26NBYvXtz5c1tbWzQ2NuZVIlSME+nH0MMBFEJuQWTatGmxY8eOLsdeeOGFOOuss7r9TG1tbdTW1uZVEgCkNWBc345XgNyCyFe+8pW4+OKLY9myZfH5z38+Nm/eHGvWrIk1a9bkdUngfWxKBsWlUpfo9iS3IHLhhRfGAw88EEuXLo1vfetbMWbMmFixYkXMnj07r0sC72NTMqDY5bqz6qc//en49Kc/neclAIASltvyXQCADyKIAADJCCIAQDK59ogAadmUDCh2ggiUMUt0gWJnagYASMaICNArNkcD8iCIAL1iczQgD6ZmAIBkBBEAIBlBBABIRhABAJIRRACAZKyaAXrFLq1AHgQRoFfsFQLkwdQMAJCMIAIAJCOIAADJCCIAQDKaVQEPtAOSEUQAD7QDkjE1AwAkI4gAAMkIIgBAMoIIAJCMZlXAc2SAZAQRwBJdIBlTMwBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwt3oGkmq78frz2fMsxx0dPHGnreagAggiQ1GvPt8TOZ15JXQaQiCACJcYIAlBOBBEoMUYQgHKiWRUASMaICNDtdE+EKR8gX4IIkHS6Z/TEkX06DpQXQQRIymgLVDZBBIrQ8aZKGieMjK+tX2gEASgrgggUoeNNlfzxwB8jwggCUF6smoES0bJzd9w863upywDoV0ZEoITsa36zV+/r66ZnPU3rmPIB8iSIQBnq6yoY0z1AKgWbmvnud78bVVVVsWjRokJdEgAocgUZEXnqqafizjvvjI9+9KOFuByUvJNdGdM44fjvG9542gnXBJCH3IPIH/7wh5g9e3bcdddd8Z3vfCfvy0FZONmpkq+tN9UClIbcg8iCBQti5syZMWPGDEEEcmardqDU5BpE7r333ti6dWs89dRTvXp/e3t7tLe3d/7c1taWV2lQNpqu/H5k2Z9GQYrxybx9XcEDVJbcgkhzc3MsXLgwNm7cGIMGDerVZ5qamuLWW2/NqyQoS92NgBSLYgxHQPHIbdXMli1bYu/evXHBBRfEgAEDYsCAAfHYY4/FD37wgxgwYEAcOXLkmM8sXbo0WltbO1/Nzc15lQcAFIHcRkQuvfTS+O1vf9vl2Ny5c2PChAlx4403Rk1NzTGfqa2tjdra2rxKAgCKTG5BpK6uLs4777wux0455ZT48Ic/fMxxAKAy2VkV6LOTaUAdOf6MGFw3uNu9ToDKUtAg8uijjxbyclARRk8cGVn253/u6X39pS8NqO+97vDG0+Jb/3FDv9UBlD4jIlDi3jsCkedy2PcuE+6L49VkSS9wlCAC9Ep/LhO2pBc4qmAPvQMAeD9BBABIxtQM0Gcn+3RggKMEEaDX/njgjxGRb1MsUFkEEaBXRk8cGa893xKrFq+L+XfMiWWzvx/N249tYO3NEl0jKsBRggjQK+8fBWnefuIrX4yoAEcJIlDG7NcBFDtBBMqY/TqAYmf5LgCQjBER4IRoOAX6gyACnBA9JkB/EESA3GmaBbojiEAZK5bpE02zQHcEEShjRhuAYmfVDACQjCACACRjagbKnEZRoJgJIlDmiqFRtFiaZoHiI4gAuTPyAnRHjwgAkIwgAgAkI4gAAMnoEYEyp1EUKGaCCJQ5jaJAMTM1AwAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQjCACACRjZ1WgW01Xfj9ee77luOdGTxxp11bgpAkiQLdee74ldj7zSuoygDJmagYASEYQAQCSMTVTgr7yyEOxc//+Y46PHzYs/uWymQkqoph019ehpwMoRoJICdq5f3/8bt/e1GVQpFL3dWhwBfpCEAG6NXriyD6fSx2EgNIiiADdMnoB5E2zKgCQjCACACRjaqYEjR82rE/HqSzd9W701O8BkIogUoIs0aUn+jqAUiKIQAX48Tfuid2v7Iv/d7A99jW/2eVcfy+pPZGVNkDlEkSgAjy14TcFW1JrRAboC82qAEAygggAkEyuQaSpqSkuvPDCqKuri4aGhpg1a1bs2LEjz0sCACUk1yDy2GOPxYIFC+LXv/51bNy4Md5999345Cc/GQcPHszzsgBAici1WfXhhx/u8vPdd98dDQ0NsWXLlvjEJz6R56UBgBJQ0FUzra2tERExrJuNt9rb26O9vb3z57a2toLUBeXOklqgWFVlWZYV4kIdHR3xd3/3d/H222/HE088cdz3fPOb34xbb731mOOtra0xZMiQvEsEAPpBW1tb1NfX9+r7u2BBZP78+bFhw4Z44oknYtSoUcd9z/FGRBobGwURACghfQkiBZmaufbaa+PBBx+Mxx9/vNsQEhFRW1sbtbW1hSip7H3lkYdi5/79xxwfP2yYLeIBKBq5BpEsy+LLX/5yPPDAA/Hoo4/GmDFj8rwc77Fz//743b69qcuAktF05ffjtedbjjne31vgA13lGkQWLFgQ99xzT/znf/5n1NXVxe7duyMior6+PgYPHpznpQH65LXnWwq2DT7wZ7nuI7Jq1apobW2NSy65JM4888zO13333ZfnZQGAEpH71AwAQHc8awYASEYQAQCSKejOqhTO+G52r+3uOFS67naYtfMs5KtgG5qdiL5siAIAFIe+fH+bmgEAkhFEAIBkBBEAIBnNqkWou+fERHhWDADlRRApQp4TA0ClMDUDACQjiAAAyQgiAEAygggAkIxm1SLU0zbsR88teviheOmt/VbRAFDSBJEi1Jtg8T97dserrW/nXwwA5EgQOQmF3O/j6AhIRMTBQ4eEEADKgiByEgq538dLb9lbBIDyo1kVAEhGEAEAkjE1UyK6W0nT0wobACh2gkiJsEQXgHIkiJyE3uz3AQB0TxA5CUYpAODkaFYFAJIRRACAZAQRACAZQQQASEYQAQCSEUQAgGQs3y0D3T0FuL+fAAwA/U0QKQOFfAowAPQnUzMAQDKCCACQjCACACQjiAAAyWhWLQPdPenXE4ABKHaCSBkolyW63S1DjrAUGaBcCSIUDcuQASqPHhEAIBkjIkXMVAUA5U4QKWKmKgAod6ZmAIBkBBEAIBlTMxSNnvY9sScKQHkSRCgamm8BKo8gUsQKOUJghQ4AKQgiRayQX/5W6ACQgmZVACAZQQQASEYQAQCSyT2I/OhHP4qzzz47Bg0aFFOnTo3NmzfnfUkAoETk2qx63333xeLFi2P16tUxderUWLFiRVx22WWxY8eOaGhoyPPSuepuhUmK1SX9VYs9PABIIdcgcscdd8S8efNi7ty5ERGxevXqeOihh+LHP/5xLFmyJM9L56qYVpj0Vy2W5wKQQm5TM4cOHYotW7bEjBkz/nyx6uqYMWNGPPnkk8f9THt7e7S1tXV5AQDlK7cg8uabb8aRI0fi9NNP73L89NNPj927dx/3M01NTVFfX9/5amxszKs8AKAIFNWqmaVLl0Zra2vnq7m5OXVJAECOcusROe2006Kmpib27NnT5fiePXvijDPOOO5namtro7a2Nq+SAIAik1sQGThwYEyePDk2bdoUs2bNioiIjo6O2LRpU1x77bV5XbbidLeixUoXAEpBrqtmFi9eHHPmzIkpU6bERRddFCtWrIiDBw92rqIpVcX05W+1CwClLNcg8vd///exb9++uPnmm2P37t3xsY99LB5++OFjGlhLjS9/AOgfVVmWZamL6E5bW1vU19dHa2trDBkyJHU5AEAv9OX7u6hWzQAAlUUQAQCSEUQAgGRybVblg139y/+IN/5woPPnE31wXjE9iA8AeksQydkHBYQ3/nCgXx5aV0wP4gOA3hJEciYgAED39IgAAMkIIgBAMoIIAJCMHpHExg3t+nyaE31eTTE9/wYAeksQydkHBYQVl/fP0lpLdAEoRRUVRLpbShuR334bAgIAdK+igoiltABQXCoqiHDi7NwKQB4EEXrlvaNJZ9efGqcMHBgREVmWsioASp0gQp+cXX9q/GrOVanLAKBMCCL0ydGRkBSNvwCUn4oKIj3tqWG/jb7R+AtAf6ioIOL/0gGguNjiHQBIpqJGRDhxR6eu3r8lPQCcDEGEXjGtBUAeBBFOiMZfAPqDIMIJMUICQH/QrAoAJCOIAADJCCIAQDKCCACQjCACACQjiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQzIDUBVSSrzzyUOzcv/+458YPGxb/ctnMAlcEAGkJIgW0c//++N2+vanLAICiYWoGAEhGEAEAkhFEAIBkBBEAIBlBBABIxqqZAho/bNgJnQOAciWIFJB9QgCgK1MzAEAyuQSRV199Na666qoYM2ZMDB48OMaNGxe33HJLHDp0KI/LAQAlKpepme3bt0dHR0fceeedMX78+Ni2bVvMmzcvDh48GMuXL8/jkgBACarKsiwrxIX+6Z/+KVatWhUvv/xyrz/T1tYW9fX10draGkOGDMmxOgCgv/Tl+7tgzaqtra0x7ANWhrS3t0d7e3vnz21tbXmXBQAkVJBm1Z07d8bKlSvjmmuu6fF9TU1NUV9f3/lqbGwsRHkAQCJ9CiJLliyJqqqqHl/bt2/v8pmWlpa4/PLL44orroh58+b1+O9funRptLa2dr6am5v7fkcAQMnoU4/Ivn374n//9397fM/YsWNj4MCBERGxa9euuOSSS+Kv/uqv4u67747q6r4NwOgRAYDSk1uPyPDhw2P48OG9em9LS0tMnz49Jk+eHGvXru1zCAEAyl8uzaotLS1xySWXxFlnnRXLly+Pffv2dZ4744wz8rgkAFCCcgkiGzdujJ07d8bOnTtj1KhRXc71ZbXw0fdaPQMApePo93ZvvvMLto/IiXj99detnAGAEtXc3HzMgMT7FXUQ6ejoiF27dkVdXV1UVVV1+762trZobGyM5ubmimtqreR7j6js+3fv7t29V45Su/csy+LAgQMxYsSID+wRLeqn71ZXV39gknqvIUOGlMQvKA+VfO8RlX3/7t29Vxr3Xhr3Xl9f36v3WcoCACQjiAAAyZRFEKmtrY1bbrklamtrU5dScJV87xGVff/u3b1XGvdenvde1M2qAEB5K4sREQCgNAkiAEAygggAkIwgAgAkU7ZB5KGHHoqpU6fG4MGDY+jQoTFr1qzUJRVUe3t7fOxjH4uqqqr4zW9+k7qc3L366qtx1VVXxZgxY2Lw4MExbty4uOWWW+LQoUOpS8vFj370ozj77LNj0KBBMXXq1Ni8eXPqknLX1NQUF154YdTV1UVDQ0PMmjUrduzYkbqsJL773e9GVVVVLFq0KHUpBdHS0hJXXnllfPjDH47BgwfH+eefH08//XTqsnJ35MiRuOmmm7r8Xfv2t7/dp2e2lYKi3ln1RN1///0xb968WLZsWfzN3/xNHD58OLZt25a6rIK64YYbYsSIEfHss8+mLqUgtm/fHh0dHXHnnXfG+PHjY9u2bTFv3rw4ePBgLF++PHV5/eq+++6LxYsXx+rVq2Pq1KmxYsWKuOyyy2LHjh3R0NCQurzcPPbYY7FgwYK48MIL4/Dhw/G1r30tPvnJT8Zzzz0Xp5xySuryCuapp56KO++8Mz760Y+mLqUg3nrrrZg2bVpMnz49NmzYEMOHD48XX3wxhg4dmrq03N1+++2xatWqWLduXZx77rnx9NNPx9y5c6O+vj6uu+661OX1n6zMvPvuu9nIkSOzf/3Xf01dSjL/9V//lU2YMCH73e9+l0VE9swzz6QuKYnvfe972ZgxY1KX0e8uuuiibMGCBZ0/HzlyJBsxYkTW1NSUsKrC27t3bxYR2WOPPZa6lII5cOBAds4552QbN27M/vqv/zpbuHBh6pJyd+ONN2Yf//jHU5eRxMyZM7MvfelLXY599rOfzWbPnp2oonyU3dTM1q1bo6WlJaqrq2PSpElx5plnxqc+9amKGRHZs2dPzJs3L/7t3/4tPvShD6UuJ6nW1tYYNmxY6jL61aFDh2LLli0xY8aMzmPV1dUxY8aMePLJJxNWVnitra0REWX3O+7JggULYubMmV1+/+XuF7/4RUyZMiWuuOKKaGhoiEmTJsVdd92VuqyCuPjii2PTpk3xwgsvRETEs88+G0888UR86lOfSlxZ/yq7IPLyyy9HRMQ3v/nN+MY3vhEPPvhgDB06NC655JLYv39/4urylWVZfPGLX4x/+Id/iClTpqQuJ6mdO3fGypUr45prrkldSr96880348iRI3H66ad3OX766afH7t27E1VVeB0dHbFo0aKYNm1anHfeeanLKYh77703tm7dGk1NTalLKaiXX345Vq1aFeecc0488sgjMX/+/Ljuuuti3bp1qUvL3ZIlS+ILX/hCTJgwIf7iL/4iJk2aFIsWLYrZs2enLq1flUwQWbJkSVRVVfX4OtonEBHx9a9/PT73uc/F5MmTY+3atVFVVRX//u//nvguTkxv733lypVx4MCBWLp0aeqS+01v7/29Wlpa4vLLL48rrrgi5s2bl6hy8rRgwYLYtm1b3HvvvalLKYjm5uZYuHBhrF+/PgYNGpS6nILq6OiICy64IJYtWxaTJk2Kq6++OubNmxerV69OXVrufvazn8X69evjnnvuia1bt8a6deti+fLlZRfCSqZZ9frrr48vfvGLPb5n7Nix8cYbb0RExEc+8pHO47W1tTF27Nh47bXX8iwxN72991/96lfx5JNPHvMsgilTpsTs2bNL8j/e3t77Ubt27Yrp06fHxRdfHGvWrMm5usI77bTToqamJvbs2dPl+J49e+KMM85IVFVhXXvttfHggw/G448/HqNGjUpdTkFs2bIl9u7dGxdccEHnsSNHjsTjjz8eP/zhD6O9vT1qamoSVpifM888s8vf84iIiRMnxv3335+oosL56le/2jkqEhFx/vnnx+9///toamqKOXPmJK6u/5RMEBk+fHgMHz78A983efLkqK2tjR07dsTHP/7xiIh4991349VXX42zzjor7zJz0dt7/8EPfhDf+c53On/etWtXXHbZZXHffffF1KlT8ywxN72994g/jYRMnz69cxSsurpkBvx6beDAgTF58uTYtGlT55L0jo6O2LRpU1x77bVpi8tZlmXx5S9/OR544IF49NFHY8yYMalLKphLL700fvvb33Y5Nnfu3JgwYULceOONZRtCIiKmTZt2zDLtF154oWT/nvfFO++8c8zfsZqams6R/7KRuls2DwsXLsxGjhyZPfLII9n27duzq666KmtoaMj279+furSCeuWVVypm1czrr7+ejR8/Prv00kuz119/PXvjjTc6X+Xm3nvvzWpra7O77747e+6557Krr746O/XUU7Pdu3enLi1X8+fPz+rr67NHH320y+/3nXfeSV1aEpWyambz5s3ZgAEDsttuuy178cUXs/Xr12cf+tCHsp/85CepS8vdnDlzspEjR2YPPvhg9sorr2Q///nPs9NOOy274YYbUpfWr8oyiBw6dCi7/vrrs4aGhqyuri6bMWNGtm3bttRlFVwlBZG1a9dmEXHcVzlauXJlNnr06GzgwIHZRRddlP36179OXVLuuvv9rl27NnVpSVRKEMmyLPvlL3+ZnXfeeVltbW02YcKEbM2aNalLKoi2trZs4cKF2ejRo7NBgwZlY8eOzb7+9a9n7e3tqUvrV1VZVmZbtAEAJaP8JtEBgJIhiAAAyQgiAEAygggAkIwgAgAkI4gAAMkIIgBAMoIIAJCMIAIAJCOIAADJCCIAQDKCCACQzP8H56ntsE6Qs5IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# mini batch算法版本\n",
    "THETA = np.random.normal(0, 0.1, size=(3,3)).reshape(3, 3) # learnable parameters #和上述公式中表述的一样，为列向量组合而成的矩阵\n",
    "EPOCH = 100\n",
    "LR = 0.5\n",
    "\n",
    "def hypothesis_batch(X, THETA):\n",
    "    X = np.reshape(X, (-1, 3))\n",
    "    #temp = np.matmul(THETA.T, x)\n",
    "    #temp = np.exp(temp)\n",
    "    # N*3 3*3  -》 N*3（N为样本的数目），3为类别数目\n",
    "    temp = X @ THETA\n",
    "    temp = np.exp(temp)\n",
    "    denominator = np.sum(temp, axis = 1)\n",
    "    hypothesis = temp / denominator[:, np.newaxis] # normalize into 1return \n",
    "    return hypothesis\n",
    "\n",
    "def predict(x, THETA):\n",
    "    x = np.reshape(x, (-1, 3))\n",
    "    h_x = hypothesis_batch(x, THETA)\n",
    "    return h_x[0]\n",
    "# def compute_loss_batch(X, Y,THETA):#X为n*3\n",
    "#     X= np.reshape(X, (-1, 3))\n",
    "#     #Y = np.reshape(y, (-1, 1))     # \n",
    "#     h_X = hypothesis_batch(X, THETA)    # hypothesis N*3\n",
    "#     #label = np.argmax(y, axis=0)  # the category of prediction\n",
    "#     loss = (-np.log(h_X[:, Y] + 0.0000001)).mean()  # loss = - y * log(y')\n",
    "#     return loss\n",
    "\n",
    "def update_parameters_batch(THETA, x, y):\n",
    "    x = np.reshape(x, (-1, 3))\n",
    "    y = np.reshape(y, (-1))\n",
    "    \n",
    "    h_x = hypothesis_batch(x, THETA)\n",
    "    # for row, index in zip(h_x, y):\n",
    "    #     loss = -np.mean(np.log(h_x[:, y]+ 0.0000001)) # loss = - y * log(y')\n",
    "    #loss = -np.log(h_x[:, y]+ 0.0000001)# loss = - y * log(y')\n",
    "    loss = np.sum(-np.log(h_x[np.arange(h_x.shape[0]), y]+0.0000001))\n",
    "    #print(f'h_x[:, y]: {h_x[:, y]}, loss: {loss} , y:{y}')\n",
    "    for k in range(3):\n",
    "        indicator = y == k\n",
    "        indicator = np.array([int(b) for b in indicator])\n",
    "        grad = (indicator-h_x[:, k])[:,np.newaxis] * x\n",
    "        grad = np.mean(grad, axis= 0)\n",
    "        THETA[:, k ] = THETA[:, k]+LR*grad\n",
    "\n",
    "    return THETA, loss\n",
    "\n",
    "for epoch in range(EPOCH):\n",
    "    #LR = LR * (1 / (1 + DECAY_RATE * epoch))\n",
    "    i = 0 # retrieve H_x\n",
    "    batch_size = 8\n",
    "    iterations = len(X_train) // batch_size\n",
    "    idxs = np.array(list(range(len(X_train))), np.int32)\n",
    "    np.random.shuffle(idxs)\n",
    "    loss = 0\n",
    "    #print(X_train.shape, Y_train.shape)\n",
    "    for iter in range(iterations):\n",
    "        X = X_train[idxs[iter*batch_size:(iter+1)*batch_size],:]\n",
    "        Y = Y_train[idxs[iter*batch_size:(iter+1)*batch_size]]\n",
    "        #print(X.shape, Y.shape)\n",
    "        THETA, loss_batch = update_parameters_batch(THETA, X, Y)\n",
    "        loss += loss_batch\n",
    "        #print(loss_batch)\n",
    "    print('[{0}/{1}] loss is: {2}'.format(epoch+1, EPOCH, loss))\n",
    "\n",
    "\n",
    "i = 0\n",
    "print('THETA', THETA, Y_test.shape)\n",
    "H_test = np.zeros((X_test.shape[0], X_test.shape[1]))\n",
    "#H_test = np.zeros([Y_test.shape[0], 1], dtype=Y_test.dtype)\n",
    "for x, y in zip(X_test, Y_test):\n",
    "    H_test[i] = predict(x, THETA)\n",
    "    i+=1\n",
    "plt.figure(1)\n",
    "x = np.linspace(-7, 4, 50)\n",
    "plt.scatter(X_test[:, 1], X_test[:, 2], c=np.argmax(H_test, axis=1), edgecolors='white', marker='s')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References:\n",
    "- http://ufldl.stanford.edu/tutorial/supervised/SoftmaxRegression/\n",
    "- "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
